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Coherence  Eandwidth  and  Pulse  Distortion 
through  Naturally  Occurring  and 
Artificially  Modified  Ionosphere 

1.  Introduction 

A  contract  with  the  above  title  was  awarded  to  the  University  of 
Illinois  to  commence  the  investigations  on  July  16,  1978.  This  contract 
has  since  been  extended  twice  to  July  15,  1982.  Therefore,  in  this 
final  report  we  will  summarize  our  work  for  a  period  of  four  years. 

This  project  is  concerned  with  the  study  of  propagation  of  radio 
signals  through  a  perturbed  ionosphere.  These  ionospheric  perturbations 
can  be  naturally  occurring,  probably  caused  by  some  plasma  instability 
mechanisms;  they  can  also  be  man-made,  for  example,  by  high  power 
transmitters,  chemical  releases  or  nuclear  detonations  in  the  atmos¬ 
phere.  We  are  not  concerned  in  this  project  with  the  exact  physical 
mechanisms  that  produce  these  perturbations  although  investigations 
into  these  mechanisms  are  Interesting  research  problems  by  themselves. 
What  is  concerned  in  this  project  are  effects  these  ionospheric  pertur¬ 
bations  may  have  on  radio  signals  propagating  through  them.  In  general, 
these  perturbations  are  highly  complicated  and  imprecisely  known.  We 
must  then  use  whatever  information  that  is  available  to  model  these 
perturbations  either  deterministically  or  stochastically.  Such  a 
modeling  effort  will  depend  on  the  interplay  of  experimental  observations 
and  theoretical  suggestions  and  deductions.  With  accumulation  of  experi¬ 
mental  data  the  models  can  be  made  more  accurate.  Of  course  better 
ionospheric  models  will  enable  us  to  simulate  more  closely  the  true 
propagation  conditions. 


2 .  Summary  of  Problems  Investigated  and  Results 

The  problems  investigated  under  the  support  of  this  project  can 
be  broadly  classified  into  three  areas:  propagation  effects  through 
ionospheric  bubbles,  propagation  effects  through  ionospheric  irregu¬ 
larities  and  experimental  observations  of  scintillation  effects. 

These  are  described  in  the  following. 

2.1  Propagation  Effects  through  Ionospheric  Bubbles 

One  of  the  most  intensely  perturbed  ionospheric  regions  by  natural 
events  is  the  region  above  the  magnetic  equator.  Various  observations 
have  shown  that  large  depletions  in  ionization,  often  called  bubbles, 
are  possible.  Inside  these  bubbles,  the  ionization  structure  is 
extremely  complex  and  have  been  variously  described  as  ionization 
walls  and  protruding  fingers.  As  far  as  their  effects  to  radio  propaga¬ 
tion  is  concerned,  the  important  realization  is  the  fact  that  the 
ionization  gradients  can  be  very  steep,  much  steeper  than  expected 
based  on  the  usual  random  medium  hypothesis.  Because  of  this  we 
cannot  use  the  random  medium  approach  to  wave  propagation  studies; 
instead  we  must  first  model  this  bubble  medium  deterministically  and 
then  study  the  propagation  effects  through  such  a  deterministic  medium. 

The  complexity  of  the  bubble  structure  renders  any  analytical  approach 
hopeless  and  the  only  recourse  is  then  to  use  a  numerical  approach. 

Several  problems  have  been  investigated  using  this  approach.  We  summarize 
our  results  briefly  below. 

To  start  the  investigation  through  an  ionospheric  bubble  we  need 
to  do  two  tilings:  (1)  A  bull  e  model  and  (2)  A  numerical  scheme  in 
solving  the  wave  equation,  ihese  two  things  are  done  in  a  report 


(see  section  4,  publication  7).  The  bubble  model  is  constructed  using 
the  in-situ  data  and  numerical  calculations  based  on  plasma  instabilities. 
The  wave  equation  is  converted  into  a  parabolic  equation  under  forward 
scatter  approximation,  which  is  then  solved  numerically  by  an  implicit 
method  generalized  from  the  Crank-Nicholson  scheme.  Some  preliminary 
results  are  also  given  in  this  report.  These  results  obtained  from 
additional  computations  are  presented  in  a  paper  (see  section  4,  publica¬ 
tion  3).  Among  other  effects,  these  numerical  computations  show  that 
the  sharp  wedgelike  structures  associated  with  bubbles  may  be  responsible 
to  scintillation  outbursts  at  GHz  frequencies  that  are  not  predicted  by 
any  statistical  scintillation  theory.  Apparently,  sharp  gradients  in 
electron  density  variations  affect  the  field  distributions  in  ways 
quite  different  from  the  way  in  which  uniformly  disturbed  irregularities 
do.  The  effects  seem  to  be  most  distinct  at  GHz  frequencies. 

Nuclear  explosions  in  the  atmosphere  are  known  to  create  violent 
perturbations  in  the  ionosphere.  It  is  thus  of  interest  to  investigate 
the  effects  such  violent  ionospheric  perturbations  may  have  on  waves 
propagating  through  them.  Based  on  information  accumulated  in  the 
open  literature  the  gross  ionospheric  features  after  the  nuclear 
explosion  are  known,  but  for  wave  propagation  studies  we  need  to  have 
fine  structures  down  to  scales  of  ten  meters  or  less.  Since  such 
fine  details  are  unavailable  it  was  suggested  to  us  to  use  the  equatorial 
bubble  data  as  a  basis  for  construction  of  a  model.  The  horizontal 
structure  is  unchanged;  it  is  still  based  on  the  in  situ  data.  How¬ 
ever,  the  important  modifications  have  been  made  in  two  aspects: 
the  vertical  extent  of  the  bubble  and  the  vertical  correlation  distance 
of  the  spiky  structures.  To  accomplish  both  of  these  aspects  we  have 
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effectively  superposed  three  bubbles  on  top  of  each  other.  In  this 
model,  the  background  ionosphere  has  a  maximum  electron  density  of 
2.8*10^“  electrons/m2  and  an  electron  content  value  of  8.9*10^  elec¬ 
trons/nr,  corresponding  to  a  slab  thickness  of  320  km.  Superposed  ja 
this  background  ionosphere  is  ix  region  of  depleted  electron  densities 
with  maximum  depletion  equal  to  6.7*10  electrons/m2,  Inside  the 
depleted  region  there  exist  sharp  gradients  that  have  vertical  corre¬ 
lation  distances  equal  to  50  km  to  100  km.  Numerical  computations  show 
amplitude  scintillations  to  be  appreciable  at  4  GHz  (S^-0.26)  and 
clearly  discernible  even  at  a  frequency  as  high  as  15  GHz  (S^e,0.0025) . 

The  dominating  component  of  the  phase  fluctuation  has  its  origin  in  the 

variation  of  the  optical  path  which  has  the  inverse  frequency  dependence, 

-2 

The  phase  departure  from  the  optical  path  has  a  rough  f  dependence. 

At  4  GHz,  the  standard  deviation  of  phase  departure  is  0.09n  radians 
or  16  degrees.  These  results  and  additional  details  have  been  published 
(see  section  4,  publication  6). 

In  a  separate  investigation  r?e  have  also  carried  out  computations 
that  indicate  how  a  pulse  would  distort  while  propagating  through  this 
bubble  medium.  To  do  this,  the  pulse  is  Fourier  decomposed  into  discrete 
frequency  components.  Tnc  parabolic  equation  is  then  solved  for  each 
component,  which  recombines  at  the  receiver  location  to  form  a  pulse. 

The  results  are  described  in  a  thesis  (see  section  4,  publication  10). 
Since  it  is  attached  as  Appendix  B,  we  shall  not  repeat  it  here. 

2 . 2  Propagation  Effects  through  Ionospheric  Irregularities 

Under  the  support  of  this  project  we  have  also  worked  on  several 


problems  using  the  random  medium  approach.  These  are  described  in 
the  following. 


When  a  pulse  propagates  through  the  random  medium,  the  pulse  is 
expected  to  be  broadened.  In  funeral  two  distinct  mechanisms  may  be 


responsible  for  the  broadening  of  a  temporal  pulse.  In  the  first 
mechanism  the  pulse  arriving  at  the  receiver  essentially  perserves 
its  shape,  but  the  arrival  time  fluctuates  from  realization  to  realiza¬ 
tion.  This  fluctuation  in  arrival  time  may  be  caused  by  pulse  wandering. 
One  can  envision  the  presence  of  large  irregularities  which  may  be 
responsible  for  the  pulse  to  travel  along  a  different  path  and  hence 
wander  as  a  function  of  time  as  the  irregularities  move  through  the 
propagation  path.  In  the  second  mechanism  the  pulse  may  be  broadened 
due  to  scattering,  especially  multiple  scattering.  In  this  case  each 
arriving  pulse  is  spread,  resulting  in  an  broadened  average  pulse. 

To  distinguish  these  two  pulse  broadening  processes:  wandering  and 
spreading,  one  needs  to  work  with  four-frequency  mutual  coherence 
function,  T^.  The  equation  for  is  derived  and  solved  for  two 
special  cases  in  a  paper  (see  section  4,  publication  1).  The  results 
show  that  pulse  wandering  is  the  important  mechanism  under  weak 
scattering  conditions  and  pulse  spreading  is  the  important  mechanism 
under  the  strong  scattering  conditions. 

The  arrival  time  of  a  radio  pulse  as  observed  by  a  fixed  observer 
can  be  defined  by  the  first  temporal  moment  or  the  "time  centroid"  of 
the  pulse.  When  this  concept  is  applied  to  wave  propagation  in  a  random 
medium,  the  pulse  arrival  time  becomes  a  random  variable  and  will  fluct¬ 
uate  about  some  mean  value,  giving  rise  to  the  phenomenon  commonly 
called  pulse  jitters.  It  is  then  of  interest  to  investigate  the 
statistical  properties  of  this  pulse  jittex,  which  can  be  done  by 
using  the  multifrequency  mutual  coherence  functions.  This  problem  is 
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the  subject  of  a  paper  (see  section  4,  publication  2)  where,  among 
other  things,  the  statistical  moments  for  the  pulse  arrival  time  are 
derived  in  the  fully  saturated  regime.  Under  this  extreme  limit,  the 
mean  arrival  time  can  exceed  the  free  space  time  delay  by  an  appre-.  ’  t  ■■■ 
value.  At  the  same  time,  calculations  indicate  that  fluctuations  about 
the  mean  arrival  time  are  very  small  when  compared  with  the  excess  time 
delay.  This  suggests  that  the  excess  delay  time  is  caused  mainly  by 
pulse  spreading  and  not  by  pulse  wandering,  in  agreement  with  calcula¬ 
tions  carried  out  in  publication  1  (see  section  4) . 

In  applications  it  has  been  found  convenient  to  describe  the  pulse 
by  its  beginning  several  temporal  moments.  Using  such  an  approach, 
the  zeroth  moment  is  related  to  the  total  energy  in  the  pulse,  the 
first  moment  is  related  to  the  mean  arrival  time,  the  second  moment 
is  related  to  the  mean  square  pulse  width,  the  third  moment  is  related 
to  the  skewness  of  the  pulse  and  the  fourth  moment  is  related  to  the 
kurtosis  of  the  pulse.  For  propagation  in  an  ionosphere  expressions 
for  these  five  moments  have  been  derived  (see  section  4,  publications 
4  and  9).  In  addition  to  the  physical  significance  of  these  temporal 
moments,  they  can  also  be  applied  to  digital  communications.  An  ideal 
pulse  in  such  a  case  may  become  distorted  and  stretched  owing  to  propa¬ 
gation  effects.  When  the  pulsewidth  is  stretched  to  occupy  an  interval 
longer  than  one  communication  bit  inter-symbol  interference  is  expected. 
Using  the  temporal  moments  an  upper  bound  of  the  energy  content  outside 
of  one  chip  can  be  estimated.  Readers  interested  in  details  should 
consult  the  original  publications,  one  of  which  is  attached  as  Appendix  A. 

2 . 3  Experimental  Observations  of  Scint illat i on  Effects 

In  addition  to  theoretical  investigations  of  the  scintillation 
problems  we  have  also  collected  some  experimental  data  for  evaluation. 


Earlier  multifrequency  scintillation  data  received  at  various  stations 
around  the  globe  by  SRI  International  have  been  analyzed  and  the  results 
have  been  reported  in  our  Scientific  Report  No.  2  (see  section  4,  publica¬ 
tion  8).  These  early  results  will  not  be  repeated  here.  More  lately  a 
joint  effort  was  wade  in  collecting  scintillation  data  at  Ascension 
Island.  These  more  recent  data  have  been  partially  analyzed  and  some 
of  the  results  are  reported  here. 

In  the  course  of  analyzing  the  GHz  scintillation  data  from  Ascension 
Island,  it  became  clear  to  us  that  there  exist  time  shifts  up  to  one 
second  or  so  between  similar  fades  of  the  L-  and  C-band  signals. 

Initial  efforts  were  concentrated  on  eliminating  possible  causes  of 
equipment  origin.  Several  possibilities  were  considered  and  finally 
with  the  test  tape  provided  to  us  by  AFGL,  it  was  concluded  that  the 
shifts  are  real  physical  phenomenon  rather  than  the  artifact  of  the 
data  recording  process.  The  data  have  been  analyzed  and  a  physical 
model  has  been  proposed  to  interprete  the  phenomenon.  The  results 
have  been  written  up  as  a  paper  to  be  published  in  Geophysical  Research 
Letters  (see  section  4,  publication  6) .  It  is  also  included  as  Appendix 
C. 

The  Ascension  Island  data  are  still  being  analyzed.  We  report 
their  preliminary  findings  in  the  following. 

Statistics  of  Multi-Frequency  Scintillation  Signals 

The  UHF,  L-  and  C-band  data  from  Marisat  collected  at  Ascension 
Island  provide  us  with  multi-frequency  data  to  study  the  statistical 
behavior  of  the  signals.  Segments  of  these  multi-frequency  data  were 
selected  for  detailed  analysis.  Fig.  1  shows  such  an  example  where 
scintillation  indices  (S^)  for  every  minute  are  presented  for  L-  and 
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C-band  signals  during  a  ten  hour  period  on  Jan.  27,  1981.  Five  segments 
of  data  with  varying  degrees  of  scintillation  activities  were  chosen 
to  study  the  power  spectra,  the  frequency  dependence  of  S^,  and  the 
coherence  time  of  the  scintillating  signals. 

Multifrequency  Power  Spectra  of  Intensity  Scintillation 

In  Figs.  2-8,  spectra  for  UHF,  L-band  and  C-band  scintillation 
signals  are  presented  for  seven  cases.  In  all  cases,  the  UHF  scin¬ 
tillations  were  saturated  while  the  C-band  scintillations  varied  from 
very  weak  to  moderate.  The  L-band  scintillations  varied  from  strong 
to  saturation.  The  spectrum  for  weak  C-band  scintillation  (Fig.  2) 
follows  the  prediction  of  weak  scintillation  theory:  a  flat  low  fre¬ 
quency  part  followed  by  a  roll-off  at  the  fresnel  frequency  with  a 
high  frequency  asymptote  sloping  at  v«“p-l,  where  p  is  the  power  index 
of  the  three  dimensional  power-law  spectrum  of  the  form  k  P  for  the 
irregularities.  The  spectra  for  strong  scintillations  at  UHF  or 
L-band  are  quite  different  from  that  for  the  weak  scintillation  case. 

In  general,  the  spectra  are  broadened  with  the  roll-off  frequency 
pushed  further  to  the  high  frequency  end.  There  is  also  evidence 
that  the  low  frequency  part  of  the  spectrum  is  enhanced.  Table  1 
summarizes  the  important  features  of  the  spectra  for  the  seven  cases. 

It  is  interesting  to  note  that  the  slopes  of  the  high  frequency 
asymptotes  are  consistently  higher  than  most  of  the  previously  measured 
values . 

Frequency  Dependence  of 

If  we  assume  the  frequency  dependence  of  the  scintillation  index 
to  be  of  the  Form 
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then  the  spectral  index  n  can  be  computed  from  the  data.  Between 
L-  and  C-bands  signal, 

s4(c) 

n  =  -2.45  log  — -  (2) 

s4(l) 

Fig;>.  9  and  10  show  two  cases  of  the  spectral  index  n  plotted  as 
a  function  of  the  scintillation  index  of  C-band  signal.  We  note  that 
n  remains  approximately  constant  for  S4(C)<0.2,  predicted  by  the  weak 
scintillation  theory.  As  the  scintillation  gets  stronger,  n  decreases 
due  to  the  saturation  effects  on  the  L-band  signal.  This  manifestation 
of  the  multiple  scattering  effect  is  more  apparent  when  one  computes 
the  spectral  index  between  UHF  and  C-band  signals. 

According  to  the  weak  scintillation  theory,  the  spectral  index  n 
is  related  to  the  power  index  p  of  the  irregularity  spectrum  by 


n 


(3) 


Therefore,  it  is  possible  to  deduce  the  power  index  p  from  the 
scintillation  indices  of  L-  and  C-bands, 

s4(c) 

p  =  4r.-2  *>  -9.8  log  -  -  2  (4) 

s4(l) 

Figure  11  shows  a  plot  of  spectral  index  n  as  a  function  of  time 
for  the  event  on  the  27th  of  January  (Figure  1).  From  the  six  hours 
of  data,  we  selecto  .  those  periods  when  S4(C)  was  between  0.05  and  0.25 
for  computing  n  such  that  toe  weak  scintillation  theory  may  be  valid. 

We  note  the  definite  trend  of  increasing  value  of  n  after  midnight 
indicating  the  steepening  of  the  irregular if;  power  spectrum.  This 
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may  correspond  to  th’.  dissipation  of  the  small  size  irregularities 
after  midnight.  The  values  of  n  and  hence  p  from  these  computations 
may  be  somewhat  lower  than  the  real  values  for  the  irregularities 
since  some  saturation  effects  may  already  be  important  in  the  L-band 
signal  even  though  only  weak  C-band  signals  have  been  chosen  in  the 
computation. 

Since  the  power  index  p  of  the  irregularity  spectrum  can  be  obtained 
from  the  slope  v  of  the  high  frequency  asymptote  of  the  intensity  power 
spectrum  as  well  as  the  spectral  index  n  for  the  frequency  dependence 
of  S^,  it  is  interesting  to  compare  the  estimated  values  of  p  from 
both  methods.  Table  II  shows  such  a  comparison. 


Table  II 


Data  Set 

s4(c)/s4(l) 

n 

P 

v(L) 

v(C) 

P 

Jan.  26 
2135-2240 

.15/. 70 

1.64 

4.56 

4.34 

4.8 

5.34- 

Jan.  26 
2310-2400 

.09/. 47 

1.76 

5.04 

5.0 

4.9 

6-5.9 

Jan.  27 
2125-2152 

.15/. 67 

1.59 

4.36 

4.2 

4.1 

5.2-5 

Jan.  30 
2213-2258 

.17/. 85 

1.71 

4.84 

5.6 

5.4 

6.6-6 

The  values  used  in  the  Table  are  averaged  values  for  the  whole 
period  of  each  data  set.  We  note  that  the  p  values  estimated  from 
the  high  frequency  slopes  of  the  ewer  spectra  are  consistently  higher 
than  the  p  values  obtained  from  the  spectral  index  n.  This  as  explained 
above,  may  be  due  to  the  fact  that  L-band  scintillations  were  too  strong 
for  weak  scintillation  theory  to  be  valid. 


Coherence  time  of  the  UHF  signal  and  the  strength  of  scintillation 


In  strong  scintillation  theory  when  the  multiple  scattering  effects 
are  important,  it  is  known  that  the  signal  becomes  decorrelated ,  resulting 
in  decreasing  coherence  time.  Indeed,  when  the  scintillation  is  ... 
into  the  saturation  regime,  the  index  is  no  longer  an  accurate  measure 
of  the  strength  of  scintillation.  On  the  other  hand,  the  signal  coherence 
time  becomes  a  good  indicator  of  how  strong  the  irregular  structures  are. 
,'igure  12  demonstrates  this  point  as  we  compare  the  coherence  time  tc 
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Conclusion 


We  have  investigated  three  aspects  of  the  scintillation  problem 
in  this  project.  They  are:  propagation  effects  through  ionospheric 
bubbles,  propagation  effects  through  ionospheric  irregularities  and 
experimental  observations  of  scintillation  effects.  In  the  first  aspect 
the  propagation  medium  is  simulated  by  using  in-situ  data  obtained  while 
traversing  an  ionospheric  bubble.  These  bubbles  have  been  observed  to 
huve  highly  complex  structures  with  very  steep  gradients.  Because  of 
the  peculiar  nature  of  the  bubble  medium  a  deterministic  approach 
rather  a  stochastic  approach  to  solving  the  wave  equation  is  used. 

The  results  show  scintillation  bursts  especially  at  GHz  frequencies 
when  the  ray  paths  intersect  the  sharp  gradient.  These  calculations 
are  further  extended  to  1  simulated  ionosphere  under  nuclear  perturbations. 
Severe  scintillations  at  frequencies  as  high  as  15  GHz  are  possible. 

In  the  second  aspect  we  are  mainly  concerned  with  the  temporal 
behavior  of  radio  signals.  In  this  vein  we  have  investigated  the 
statistics  of  the  pulse  arrival  time  (commonly  called  pulse  Jitters), 
the  pulsewidth  and  their  effects  to  digital  communication.  Even  though 
our  study  is  theoretical  in  nature,  its  results  have  applications  to 
satellite  based  communication  and  navigation. 

The  third  aspect  of  this  project  is  concerned  with  investigations 
of  scintillation  effects  on  transionospheric.  experimental  data.  The 
experimental  data  re  recordings  made  at  various  ground  stations  of 
radio  transmissions  from  the  orbiting  DNA  Wideband  and  also  recordings 
made  at  Ascension  Island  during  a  special  campaign.  Both  sets  of  data 
are  of  high  quality,  digitized  on  magnetic  tape  for  computer  processing. 
These  data  possess  multi-channels  on  several  frequencies  and  phase 


Publications 


4.1  Journal  and  Symposium  Articles 

1.  C.  H.  Liu  and  K.  C.  Yeh,  Pulse  spreading  and  wandering  in  random 
media.  Radio  Sci.,  14(5) ,  925-931,  1979. 

2.  C.  H.  Liu  and  K.  C.  Yeh,  Statistics  of  pulse  arrival  time  in 
turbulent  media.  J.  Opt.  Soc.  Am.,  70(2) ,  168-172,  1980. 

3.  A.  W.  Wernik,  C.  H.  Liu  and  K.  C.  Yeh,  Model  computations  of  radio 
wave  scintillation  caused  by  equatorial  ionospheric  bubbles. 

Radio  Sci.,  15(3),  559-572,  1980. 

4.  C.  C.  Yang  and  K.  C.  Yeh,  Temporal  behavior  of  pulses  after 
propagating  through  a  turbulent  ionosphere.  In  "Effect  of  the 
ionosphere  on  radiowave  systems"  edited  by  John  M.  Goodman, 

pp. 525- 333,  U.S.  Government  Printing  Office,  Washington,  D.C.,  1981 

5.  K.  C.  Yeh  and  C.  H.  Liu,  Simulated  propagation  effects  on  trans- 
ionjspheric  radio  waves.  In  "Effect  of  the  ionosphere  on  radio¬ 
wave  systems"  edited  by  John  M.  Goodman,  pp. 591-598,  U.S.  Govern¬ 
ment  Printing  Office,  Washington,  D.C.,  1981. 

6.  S.  Franke,  J.  Austen,  A.  W.  Wernik  and  C.  H.  Liu,  Systematic 
refraction  caused  by  equatorial  plasma  bubbles  observed  in  micro- 
wave  scintillations.  To  appear  in  Geophys.  Res.  Letters,  1982. 

4.2  Scientific  Reports 

7.  A.  W.  Wernik,  Model  computations  of  radio  wave  scintillation  caused 
by  equatorial  bubbles.  Scientific  Report  No.  1,  AFGL-TR-79-0165 , 
July  1979. 

8.  S.  B.  Quinn,  Jr.,  Studies  of  transionospheric  scintillations  using 
orbiting  satellite  data.  Scientific  Report  No.  2,  AFGL-TR-80-0092 , 


April  1980. 


4.3  Theses 


9.  C.  C.  Yang,  Temporal  behavior  of  pulses  after  propagating  through 

a  turbulent  ionosphere.  M.S.  thesis,  1981.  Attached  as  Appendix  A. 

10.  M.  R.  Tucker,  A  deterministic  study  of  pulse  propagation  in  an 
electron  bubble  medium.  M.S.  thesis,  1982.  Attached  as  Appendix  B. 


1(, 


Figure  Captions 


Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 

Figure 


1.  Example  of  L-band  (a)  and  C-band  (b)  scintillation  data. 

2.  Intensity  power  spectrum  for  UHF,  L-band  and  C-Land  signa: 
Case  I. 

3.  Intensity  power  spectra  for  UHF,  L-band  and  C-band  signal- 
Case  II. 

4.  Intensity  power  spectra  for  UHF,  L-band  and  C-band  signal- 
Case  III. 

5.  Intensity  power  spectra  for  UHF,  L-band  and  C-band  signal- 
Case  IV. 

6.  Intensity  power  spectra  for  UHF,  L-band  and  C-band  signal- 
Case  V. 

7.  Intensity  power  spectra  for  UHF,  L-band  and  C-band  signal- 
Case  VI. 

8.  Intensity  power  spectra  for  UHF,  L-band  and  C-band  signal- 
Case  VII. 

9.  Spectral  index  n  computed  from  L-  and  C-band  data. 

10.  Spectral  index  n  computed  from  L-  and  C-band  data. 

11.  Evolution  of  the  spectral  index  n  as  a  function  of  time 

during  the  event  on  January  27,  1981. 

12.  (a)  Coherence  time  for  UHF  signal  (b)  Simultaneous  scin¬ 
tillation  index  for  C-band  signal. 
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Abstract 

A  radio  signal  after  propagating  through  a  turbulent  ionosphere 
will  suffer  distortion  owing  to  dispersion  and  random  scattering. 

For  coarse  description  of  a  temporal  signal,  the  temporal  moments 
have  been  found  to  be  convenient.  Past  studies  have  shown  that  the 
zeroch  moment  is  related  to  the  total  energy  in  the  pulse,  the  first 
moment  is  related  to  the  mean  arrival  time,  and  the  second  moment  is 
related  tv  the  mean  square  pulse  width.  In  this  report,  we  extend  the 
analysis  to  the  third  and  fourth  moments  which  are  shown  to  be  related 
to  the  skewness  and  kurtosis  of  the  pulse.  The  physical  meanings  of 
these  quantities  are  explained.  The  numerical  values  of  these  quantities 
under  typical  ionospheric  conditions  are  calculated. 

For  application  to  dig.'tal  communications,  it  is  desirable  to  know 
the  energy  content  of  the  received  pulse  outside  of  the  original  pulse 
width.  This  is  because  an  ideal  pulse,  due  to  propagation  effects  can 
be  distorted  and  stretched  to  produce  a  long  tail  occupying  an  interval 
longer  than  one  communication  bit.  Using  the  temporal  moments  the  upper 
bound  of  the  energy  content  outside  of  one  chip  can  be  estimated.  The 
numerical  values  of  these  bounds  are  also  computed. 
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1. 


Introduction 


When  signals  propagate  through  a  dispersive  and  random  medium,  they 
suffer  distortion  due  to  frequency  dispersion  and  random  scattering.  The 
main  feature  of  this  distortion  is  the  broadening  of  the  signal.  In 
binary  communication,  the  broadened  signal  in  one  bit  may  extend  into 
the  neighboring  bits.  If  it  is  of  sufficient  amplitude,  this  extension 
can  produce  serious  errors  in  the  decision  making  process.  A  practical 
example  is  the  satellite-earth  communication  link.  The  signal  in  this 
case  must  traverse  the  ionosphere  and  may  become  broadened  because  of 
dispersive  and  scattering  effects  taking  place  in  the  ionosphere.  In 
this  report,  it  is  proposed  to  describe  the  temporal  distribution  of  the 
signal  intensity  in  terms  of  its  temporal  moments  and  snow  how  these 
temporal  moments  evolve  as  the  signal  propagates  through  a  turbulent 
plasma  medium.  The  results  of  this  research  may  help  to  determine  under 
what  conditions  the  decision  making  errors  can  be  considerably  reduced. 

The  theory  of  using  first  few  temporal  moments  to  describe  the  tem¬ 
poral  distribution  of  the  received  signal  has  been  formulated  by  Yeh  and 
Liu  [1977a] .  The  works  done  before  include  the  calculation  of  up  to  the 
third  moment  in  a  nondispersive  medium  [Liu  et  al.,  1978]  and  up  to  the 
second  moment  in  a  dispersive  medium  [Yeh  and  Yang,  1977c;  Liu  and  Yeh, 
1977;  Yeh  and  Liu,  1977b;  Yeh  and  Liu,  1979].  In  Chapter  2  of  this  report 
the  evaluation  of  temporal  moments  is  extended  to  the  fourth  order.  Mean¬ 
while  the  medium  will  be  assumed  to  be  dispersive  as  well  as  turbulent. 

In  this  calculation,  the  two-frequency  mutual  coherence  function,  T ,  is 
expanded  into  a  power  series  of  the  relative  difference  of  the  wave 
numbers.  For  the  purpose  of  computing  temporal  moments  it  is  sufficient 


to  evaluate  the  series  coefficients  only  up  to  the  desired  order. 

The  equation  for  r  will  be  derived  by  using  the  Markov  approximation 
under  the  assumption  of  forward  scattering. 

Chapter  3  contains  a  discussion  of  the  physical  significance  of 
the  moments  formulated  in  Chapter  2.  Quantities  such  as  skewness  and 
kurtosis  will  then  be  defined  in  terms  of  the  temporal  moments.  These 
quantities  ere  used  to  describe  the  temporal  characteristics  of  the 
signal  intensity  distribution.  The  formulas  will  then  be  specified  by 
two  sets  of  parameters  that  describe  the  geometry  and  content  of  the 
medium.  Numerical  results  will  be  presented  as  graphs. 

Si  ae  properties  of  moments  will  be  used  in  Chapter  A  to  describe 
the  trailing  tail  of  the  received  signal.  By  using  these  properties, 
the  upper  bound  for  the  fractional  energy  of  the  signal  outside  a  time 
interval  can  be  calculated.  With  binary  communication  application  in 
mind,  we  calculate  the  upper  bound  of  the  fractional  energy  of  the  signal 
contained  in  the  neighboring  bits.  Some  numerical  results  will  be  shown. 

The  geometry  of  the  problem  is  shown  in  Figure  1.  We  assume  that 
there  exists  a  plasma  medium  with  a  homogeneous  background  in  the  region 
z>0  and  that  random  irregularities  exist  only  inside  the  slab  between 
z-O  and  z«L.  The  positions  of  the  transmitter  and  the  receiver  are  also 
shown  in  this  figure.  Note  that  the  receiver  is  always  at  a  position 
with  z^L.  Carrier  plane  waves  with  a  Guassian  envelope  are  assumed  to 
be  incident  at  z-0. 


3(> 


TO  TRANSMITTER  AT  -oo 


W777J7I777777T77777T77TT7777, 

RECEIVER  AT  (0,0, z) 


Figure  1.  The  geometry  of  the  problem.  Note  that 
plasma  occupies  the  whole  region  0<z'^_z, 
but  only  inside  the  slab  CKz'<L,  there 
exists  random  irregularities. 
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2. 


calculation  of  Temporal  Moments 


We  give  the  definitions  of  the  moments  first.  The  equation  for 
two-frequency  mutual  coherence  function,  F,  is  then  derived  under  the 
forward  scattering  assumption  and  Markov  approximation,  r  is  expanded 
into  a  power  series  of  the  relative  difference  of  wave  numbers.  The 
coefficients  of  this  series  up  to  fourth  order  are  evaluated.  Moments 
up  to  fourth  order  can  then  be  computed.  This  method  for  calculation 
of  temporal  moments  was  Introduced  by  Yeh  and  ^iu  [1977a]. 

The  computations  for  the  third  and  the  fourth  moment  are  quite 
tedious.  We  can  imagine  that  the  work  of  evaluating  the  fifth  or  any 
higher  order  moment  will  be  more  tedious  and  it  will  be  easy  for  us 
to  make  mistakes.  So  we  will  stop  at  the  fourth  moment,  although  any 
higher  order  moment  can  give  us  more  information  about  the  signal 
shape. 


2.1  Definition  of  the  Moments 

The  problem  of  plane  wave  propagation  in  a  random  medium  is  usually 
formulated  by  the  equation 


P(z,t) 


f  (cu)u(z,u))e 


j  [ut-k(cj)  z] 


d(D 


(1) 


Here  we  assume  that  the  wave  propagates  in  the  z  direction.  In  this  equa¬ 
tion,  u(z ,uj)  is  the  complex  amplitude  and  is  assumed  to  be  unity  at  the 
boundary,  i.e.,  u(0,w)“l  at  z*0  in  our  geometry  of  the  problem.  u(z,w) 
is  used  to  describe  the  random  effect  on  the  signal  component  at  circular 
frequency  w  In  the  medium.  In  general  u  is  a  random  function  of  transverse 

■f 

coordinate  p  as  well  as  z  and  w.  When  there  is  no  cause  of  confusion,  the 
explicit  indication  of  such  dependences  will  be  suppressed.  f(u.)  is  the 
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frequency  spectrum  of  the  signal  impressed.  k(w)  is  the  wave  number 
and  usually  includes  the  dispersive  characteristics  of  the  medium. 

We  note  for  real  p(z,t),  both  f(aj)  and  u(z,w)  are  required  to  be  even 
in  to,  which  we  assume. 

Rewrite  equation  (1)  in  the  form 


P(z,t)  *  ReA(z,t)  exp [ j (wct-kcz) ] 


(2) 


This  represents  a  wave  of  carrier  frequency  to  ,  wave  number  kc*k(toc) , 
and  slowly  varying  complex  envelopa  A.  A  is  given  by 


00 
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A(z,t)  - 

—00 


j [ftt-(k-k  ) z] 

F(fi)U<at,n)e  C  dfi 


(3) 


where  F(ft)*f  0oc+f2) ,  U(z,Q)*u(z,<oc+Q)  and  .  Note  that  A(z,t)  is 

random  since  u(z,w)  is  random. 

Define  the  nth  temporal  moment  by  the  equation 


M(n)(z)  = 


n 


<A*(z,t)t  A(z,t)>dt  n*0,l,2, 


•  •  •  • 


(A) 


Here  the  notation  <  >  denotes  an  ensemble  average.  Insertion  of  equation 
(3)  in  (A)  leads  to 


M(n)(z) 


rr 

j 


— OQ 


F*(n2)F(fi1)rtn  exp{ j [fi1-n2)t-(k1-k?)z] }dfi2dft2dt 


(5) 


where  ki=k(o)i)  and  k2=k(w2)  and  T  is  the  two-frequency  one-position  mutual 
coherence  function  given  by 


r  =  <u,p,z,ni)u*(o,z,n2)>  **  <u(o ,z,idi)u*(p ,z,w2) 


In  order  to  proceed  further  In  equation  (5)  we  need  the  relation 
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where  6^n\‘fli-Q2)  is  the  nth-order  Dirac  delta  function.  Using  equation 
<’6),  equation  (5)  becomes 


M(n)(z)  -  2TT(j)n 


F*(fi2)e 


jk2z  3n[F(^L)r  e"jklZ] 


dO- 


an, 


(7) 


Further  integration  of  (7)  requires  knowledge  in  T.  The  form  of  (7) 
suggests  an  expansion  of  r  in  the  series 


r<klfk2)  •  r0+r1n+r2  n2+  .... 


(8) 


where  n*(k2“ki)/k2»  and 


r  - 


i_  aT 
n:  an11 


(9) 


n-0 


As  can  be  seen  in  (7),  the  nth  temporal  moment  depends  only  on  T0,  Flf  ... 
rn  and  not  on  F^^  and  higher  order  terms. 

Once  To,  ...  r*n  arc  calculated,  we  can  obtain  M^(z)  quite 
straightforwardly  although  the  operation  can  be  very  tedious  as  we 
shall  find  later  on. 


2 . 2  Equation  for  Two-Frequency  Mutual  Coherence  Function 

To  derive  the  equation  for  F,  let  us  start  from  the  Helmholtz  wave 
equation  for  the  wave  function  41 

v2'B-k2u+eO't'  -  o  do) 

where  the  time  dependence  exp(jut)  is  understood.  The  random  function  K 
is  assumed  to  be  a  homogeneous  random  field  and  is  independent  of  frequency. 
The  frequency  dependence  is  included  in  the  factor  8.  In  a  turbulent  plasma 
medium 


40 


8  •  - 


I 


5 


AN 

N 


(11) 


where  w  is  the  circular  plasma  frequency  and  N  is  the  electron  density 
of  the  background  medium  which  is  assumed  to  be  homogeneous.  AN  is  the 
fluctuating  electron  density.  The  validity  of  equation  (10)  for  various 
waves  has  been  discussed  by  Tatarskii  [1971].  In  general  it  requires 
that  the  typical  turbulent  scale  be  large  compared  with  the  wavelength, 
which  we  assume. 

From  equation  (10)  and  under  the  forward  scattering  assumption  we 
can  obtain  the  parabolic  equation  for  the  complex  amplitude  u  as 

V  2u-2jk  +  k2S(wK(p,z)u  -  0  (12) 

1  a  Z 

,  g  2  g  2  + 

where  u  and  V  are  related  by  Y-u  exp(-Jkz)  and  +  jp?  and  P*(x,y). 

From  equation  (12),  under  the  Markov  approximation,  the  equation  for  the 
two-frequency,  two-position  mutual  coherence  function 


r(Pi,P2>  =  <u(pi  ,Z,wi)u*(p2,z,aj2)> 


can  be  derived  [Tatarskii,  1971] 


sr(oi  ,p2)  J 


dZ 


2kik2 


Ti  1  Tf 


+  T^-r-  (k2V2-kiV2  )r(p!,p2)  -  J  A  (p2-o1)r(p1,p2)  -  0  (13) 

£  P 


Here  we  denote  this  two-position  function  by  F^j.pn),  so  that  it  can  be 

distinguished  from  the  one-position  (pi_p2)  function  F  in  Section  2.1. 

The  dependence  on  frequencies  is  suppressed  when  it  is  not  explicitly 

needed.  In  the  equation  above  A  (po-pi)  is  givin  bv 

P 


Ap(P2-Pi) 


-1 


[  (ki  28i  2+k2tB22)AJ.  (0)  -  Ikji-ijk  P.  A.-(P2-f*i)  ] 


where  Bi  =  c(«*>i),  B2~B(cu2)  and  Af(p)  is  defined  by 


■1  l 


(14) 


A  ^(p)  = 


B^(p,z)dz 


(15) 


and 


B^Cp.z)  =  <5(pV  ,z'K(p '+p,z'+z)>  (16) 

where  £(p,z)  Is  a  homogeneous  random  field.  If  plane  waves  are  Impressed 
ac  z*0,  r(pltp2)  ia  a  function  of  P2-Pi"P  and  equation  (13)  becomes 


ar(p)  ,  JAk_ 

3z  2kik2 


vT2r(p)  -  A  a  (p)r(p)  -  o 


(17) 


where  Ak*k2-ki,  The  boundary  condition  for  (17)  is  F(p)-1  at  z»0. 

Since  for  the  calculation  of  up  to  the  fourth  moment,  only  Fq,  F}, 
r_,  r 3 ,  r4  are  required,  it  is  not  necessary  to  solve  equation  (17)  for 
T  itself,  but  its  expression  suggested  in  (18). 

When  0<z<L,  let's  assume 


F(p)  -  W(z)exp-Hz) 


(18) 


Where 


♦(*)  -  *(ki2Bi2  +k22B22)A^(0)z/8 


Expand  W(z)  in  the  form 


w  -  W0+wlrrfW2n2+ 


(19) 


If  equations  (18)  and  (19)  are  substituted  in  (17),  inside  the  turbulent 
slab,  i.e.,  0£Z<L,  we  obtain  the  following  general  equation: 

“il  “  4  k^B2^A£;(p)Wn  -  I"  2k7  7X‘  +  ^  k22S22A,(p))  [W^+W^*  ...  +W!+W0] 

n  ■  1,  2,  3,  .... 


4  P 


(20) 


and 

3W0  1 

57-  "  4  k22622A^(p)W0  -  0  (21) 


The  boundary  conditions  at  z*«0  are 


W0-l,  W1-W2"  _ -  0  (22) 

Note  that  in  the  process  above  we  have  used  the  expansion  of  Bi 

81  -  82^282 n+3 62 n+  ...  +  (n+1) B2nn+  .  (23) 

which  can  be  shown  to  be  valid  by  using  (11).  When  z>L,  i.e.,  outside 
the  turbulence  slab,  we  write 

T(p)  -  W'(z)  exp$(L)  (24) 

We  also  expand  W'(z)  into  a  power  series 

W'(z)  -  Wo'-Wi  'n+W2'n2+  .  (25) 


The  corresponding  equations  of  (20)  and  (21)  can  be  written  directly  by 
setting  A^(p)-0  in  equations  (20)  and  (21),  i.e,, 


3W ' 

_ r 

3z 


2k2  V(,,n-lW'.-2(  ••• 

n  ■  1,2,3,  .... 


and 


(26) 


3W0* 

—  "°  (27) 


The  boundary  conditions  at  z«L  are 


W0 1 (L)  -  W0(L),  W, '(L)-Wi(L) , 


/n 


(28) 


From  equation  (20),  we  rewrite  the  same  equation  for  i.e. 


rfW  ,  1 


Vp)Wn-l 


-J 


r -  o  <•+ 

l2k2  VT  4 


2+  -  k22a-2 


■&2' 


Vp) 


^^n-2+Wn-3+‘ 


Subtracting  (28)  from  (20)  we  obtain 


3(W  -tf  .)  1 

n  n-i 


~  4  k22S22Ac(p>  (Wn-Wn-1)  “  I  ft  V+  I  "•rer/Af.(p)]Wn_1 


Then  a  set  of  recurrent  formulas  for  V  can  be  obtained  for  0<z<L 

n  —  _ 


W  (z)  -  W  (z)  + 
n  n-i 


z  1 
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£  k22822A^.(p)(z-2') 


I-  2kT  V*  i  k22B22A?(p)] 


Un  A z')dz' 
n-i 


n«2,3,4 


z  1 


W;(z)  -  W0 (z)  + 


:  4  k22022AF(p)(z-z')  ,  ^ 

*  ’  t_  2kJ  V+  T  k22B22A^(p)]W0 

•“  k22822A  (p)z 

-e  ^ 


It  is  easy  to  obtain  a  similar  set  of  >.,-ant  formulas  for  W  * 

n 


Wn*(z)-Wn-(L)  -  IWn.1(«)-Wn.1(L)]  -  2tT  J 

L 

n»l,2,3,  ... 


YVi(z’>dz 


Note  that  in  deriving  equation  (30),  (31),  and  (32)  we  have  used  the 
boundary  conditions  in  equations  (22)  and  (28).  From  equations  (21) 
(22),  (26)  and  (27)  we  can  evaluate  Wq  and  W0 ' 

W0(z)  -  exp  [  ^  k2:^:.2Ar  (p)z] 


.+W;+Wa1 

(29; 

(30) 

(31) 

•  • 

(z')dz' 

(32) 

f 

(33) 

9 

(34) 


AA 


and 


Wo' (2)  -  W0(D  -  exp  [  ^  k.22622Aj.(D)L] 


(35) 


Using  equation  (34)  and  (35)  as  the  starting  point,  the  recurrent  formulas 
(31)  and  (33)  will  lead  to  Wr(z)  and  then  wn'(t)  any  n»l,2,3,  ...  . 
After  we  have  these  results,  equations  (23),  (24)  and  the  relation 


„  i_  ar  (p) 


n  n 


l  3n° 


n*0,  p-0 


n»0 ,1,2,  ... 


(36) 


can  be  used  to  obtain  Tq*  Fj,  f2,  Tj  and  r4  for  2>L. 


2 . 3  The  Formulas  for  Tn.  Ti.  f?,  Fq,  Fu 

If  we  assume  the  random  field  £*AN/N  is  Isotropic  as  well  as  homo¬ 
geneous,  then  A.(p)«A  (p).  Since  A  (p)  is  an  even  function,  we  can  ex- 
pand  A^(p) 

A^(p)  ■  A^(p)  -  Ao+A2 32+A4P4+Aep6+A8P8+  .....  (37) 

So 

Vp)|„-0  ■  A» 


VtVo)|p-0  '  0 


VT2Ae(o) 


3-0  '  4A2 


Having  these  identities  above,  we  are  ready  to  calculate  F q ,  F j ,  f2,  F 3 ,  F 4 . 
As  stated  in  the  last  paragraph  of  Section  2.2,  once  W0',  Vt ' ,  V.\  '  ,  W?',  W4', 
are  obtained,  we  can  evaluate  Tq,  r 1 ,  F2,  F 3 ,  r4  from  equation  (24)  and  (36). 
For  the  region  z>L,  the  results  are  as  follows: 
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(38) 


Ti  -  -  |  k2822A2(2Lz-L2)  (39) 

r2  ■  -  ^  k2‘$22AoL  -  ^k2622A2 (2Lz-L2) 

-  822A4(3Lz2-3L2z+L3)-k22824A22(i  L2z2-  L3z  +  ■—  L4)  (AO) 

1*3  -  -  ^  k22S22A0L  -  k2822A2 (2Lz-L2) 

622 

-  2022A4(3Lz2-3L2z+L3)+j  —  A6(12Lz3-18L2z2+12L3z 
-3L4)-k22624A22(L2z2-  |  L3z+  |  L4) 

+  ifeo  k23^26A23(720L3z3-1320L4z2+8A6L5z-185L6) 


+  k23024  ApA2 (2L2z-L3)4-  £  k2024A2A4(18L2z3 

-  31L3z2+20L4z  -  L5) 

r4  "  Tfa  k24e24A02L2  -  |  k22622A0L 


+  k2482SA24(Y^  L4z 4  -  ~  L5z3  +  llz2  -  -U21-  t.7«  4.  -1212-  T  8 
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IIL-l oz,  . 

5760  L  z  806A0 


L'‘  +  64tnoL*> 


A  jk23fa2bA,3(|  L3z3-  L4Z2+  L5z  -  ^  Lfc) 
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-k22024A22(|  L2z2-  y|  L3z+  y|  L4) 


- j  k2  02  2A  2 ( 2Lz-L2 ) 


+k24B26A0A22(±?  LV-  y^  L4z  +  ^  l5) 

+jk2  3824  A0A2(i  L2 z-  ±  L3) 

+824A42(12L2z‘*-28L3z  3+  |1  Luz2-  i|i  L8Z+  |0  L&) 
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4  r. 


-622A4(12Lz2-12L2z+4L3) 


-jk2024A2A4(15L2z3-  L322+  |0  l4z_  ^5) 

+k22B26A22A4(3L3z4-  f  L**z3+  L V.  f|f  L^z*  L7} 

+k22B24A0A4(|  L2z2-  -|  L3z+  f  L4) 

<*22 

+J  A6(48L«3-72L2z2+48L3*-12L4) 

+024A2A6(f  L2z4-51L3z3+  if  L4z2-  if  L5z+  f  L6) 

e2? 

+  A8(96Lz4-lQ2L2z3-t-192L3z2-96L4z+  f  Ls) 


(42) 


where 


k2  •  (tu2/c)  ( 1— ujp 2  /u,22)1/2 


(43) 


and 


02  ■  -U)  2/(oj22-u>  2) 

p  p 

in  a  turbulent  plasma. 

2 • A  Formulas  for  Temporal  Moments 

Now,  we  consider  a  plane  wave  Gaussian  pulse  at  carrier  frequency 
uc  PrcPa8atiufe  through  a  turbulent  plasma.  The  envelope  of  the  original 
pulse  is  given  by 


A(t)  -  exp(-t2/T02) 


In  other  words,  we  have  the  frequency  spectrum 


f  fa) 


2  it 


-t2/T02  jwt 

e  u  cosuj  t  eJ  dt 
c 

O 

+exp  [  -T0 (U)—A,c ) 2  /4  ]  >  4y 


—  (exp[-Tg':(o-f-ty  )2/4] 
4  ►  '  TT  C 


(44) 


From  the  relation  F(Q)«f (wc+ft)  and  we  have 

F(fl)  -  F_(fl)+F+(n) 


where 

To 

F-(fl)  = -  {exp[-T02(W-  2w  )2/4]} 

4^ 

To 

F+(fl)  2 -  {exp[-To2ft2/4] } 

4/tt" 

Since  the  value  of  Tq  gives  an  order  of  magnitude  for  the  original 
signal  width,  usually  Tq^  is  much  greater  than  one,  which  we  assume. 
This  implies  that  F-(fi)  and  for  n-1,2,3,4  are  nearly  zero 

in  the  region  fl?" -w  ,  meanwhile  F  (£1)  and  for  n*l,2,3,4 

have  negligible  values  in  the  region 

Now  let  us  divide  the  integration  in  equation  (7)  into  two  parts 

as: 

M(n)(z)  -  M_(n)(z)+M+(n)(z) 


where 


-Cd 


M_(n)(z)  =  2tr(  j  )n 


F*(n2)ejk2Z  antF(n1)re~:jklZ]/a^1n 


n  i  2 


and 


M+(n)(z)  5  2 tt ( j ) n 


F*(n2)eJk2Z  3niF(n1)rejK'‘Z  ]/3fiin| 


1  «l-«2 

dSl2 


Thun  for  the  reasons  stated  in  the  last  paragraph,  we  can  approximate  , 
M_^n^(z)  and  M_^n^(z)  as 


to 


M_(n)<z)  =  2ir(j)n  F*_(n2)eJk2Z  an[F„(«1)re"jklZ]/3a1n| 

1  'ai-fl2 

dft2 


and 


(45) 


M  (n)(z)  5  2ir(j)n  F*  (n2)eJk2Z  3°[F  <Ox)r«"^kl*] /afii” 

dn2 

(46) 

If  we  calculate  the  partial  derivative  inside  the  Integrand  In  equa¬ 
tion  (45)  and  use  the  results  in  equation  (38)  through  (42)  for  F0,  r x ,  r2,  r3, 
r4,  we  can  express  M_^(z)  as  a  summation  of  terms  which  have  the  general  form 


* 

F*.(n2)OmF.(n2)/9n2m)g(z,n2)dn2  (47) 

—40 

where  m  is  one  of  0,  1,  2,  3,  4.  In  the  last  integration,  g(z,ft2)  has  the 
form  as 

3k2  j  32k2  k  33k2  l  3uk2  m 
g(z,n2)  -  B 82*^2  (g^z-)  (3'jj2j)  (^7) 

where  B  is  some  constant  in  terms  of  L,  z,  Aq»  A2,  A4,  Ag  and  Ag,  and  h, 

i»  j»  k,  l,  m  are  some  nonnegative  integers.  Note  that  from  the  expressions 
for  B2  in  equation  (11)  and  k2  in  equation  (43),  we  can  find  that  g  ( z ,  Si  2 )  Is 

either  an  even  function  or  an  odd  function  of  w2,  since  Q2“U)2“U)C- 

Doing  the  same  thing  for  M  ^(z)  in  equation  (46),  we  can  obtain 
an  identical  summation  expression  for  M+^(z)  except  that  the  corres¬ 
ponding  general  expression  of  (47)  becomes 

QO 
' 

F*+(n2)(3mF+(n2)/3n2m]g(z,n2)dp.; 

—00 

which  just  has  different  spectrum. 
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(48) 


For  the  purpose  of  further  calculation  in  expressions  (47)  and  (48) , 


we  expand  g(z,H2)  into  a  power  series  in  two  ways  as 

g(z.^2)  ■  80+81^2+82^2 2+  •• 


and 


g(z,02)  ■  go'+gi '(^2+2w  )+g2' (^2+2“)2+  ••• 

c  c 


where 
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a«2  '^2"0 


£-0,1, 2, 3,  .... 


and 


l  3  g(z,fl2> 

^  3fl2l 


0->»-2u) 


£-0,1, 2,3,  .... 


0 

note  that  if  g(z,il2)  ia  even  with  respect  to  02--u)c(cu2"0) ,  g&— (—1)  g^ 

i+1 

and  if  g(z,il2)  is  odd  with  respect  to  fi2"_u)c»  (— 1)  ' .  Having 

these  results,  we  can  rewrite  expression  (47)  as 
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F-*(02)[3raF.(02)/302n’]g(2,n2)dn2 

— oo 


P.*(n2)  1 3mF_(«2)/^2m]  [go  ’-fgi '  (ft2+2wc)+g2  '  (fi2+2u)c)  2+  .  .  .  .  ]dn2 


F+*(fto  ' )  [  3mF^(fl2  ' )  /  3il2  ,ir]  [go '+81  '02  '+g2  'n2 ' 2  +  .  ...]dn2( 


where  we  have  replaced  the  variable  02  by  02  '-02+2ojc. 
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On  the  other  hand,  using  equation  (49a)  we  can  also  rewrite 
expression  (48)  as 

00 

F+*(fi2)0111F+(n2)/3fi2m)g(z,fi2)dfi2 

—00 

00 

f 

■  F+*(fl2) (3mF+(n2) /3n2a) [gQ+giQ2+g2n22+. . . ]dft2 

— 05 

(51) 


In  the  first  order  calculation,  we  neglect  g2'ft2'2  term  and  all  higher 
order  terms  in  equation  (50)  and  neglect  g2ft22  term  and  all  higher  order 
terms  in  equation  (51).  Now,  if  we  change  the  dummy  variable  '  into 
ft2  in  equation  (50)  and  add  it  to  equation  (51)  in  both  sides,  we  can 
conclude  that  (z)“M_^n^ (z)+K+^n^ (z)  can  be  approximated  by  a  summa¬ 
tion  of  terms  with  a  general  form 
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F+*(n2)(3mF+(n2)/3fi2m)g(z,n2)dn2 
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if  g(z,ft2)  is  even  with  respect 
to  w2-0  and  m  is  even 
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F+*(n2)[3nV+(si2)/3fi2m]n2dn2 


if 

g(z,n2) 

is 

even  with 

respect 

to 

n)2*0  but 

m 

is  odd 

if 

g(z,il2) 

is 

odd  with 

respect 

to 

u2-0  and 

m 

is  odd. 
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if  g(z,fi2)  is  odd  with  respect 
to  oj2«0  but  m  is  even 


So,  in  summary,  for  the  purpose  of  obtaining  the  formula  for  the 
nth  moment,  we  first  expand  equation  (43)  and  (46)  to  get  summations  of 
terms  with  general  forms  in  expression  (47)  and  (48)  respectively.  Then 


!>1 


equation  (52)  is  used  to  reduce  each  ten?  into  a  simple  integration 
form  like  that  gi  "jn  by  the  right  hand  side  of  equation  (52).  After 
finishing  those  simple  integrations,  we  add  them  up  to  obtain  the 
formula  for  the  nth  moment.  In  order  to  shorten  those  lengthy  for¬ 
mulas,  we  use  the  notations  X-oi  2/w  2  and  Z-l-X.  Those  formulas  are 
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In  order  to  understand  the  physical  significance  to  these  moments 
we  shift  the  origin  of  time  to  M(1)/M(0),  which  will  be  defined  to  be 
the  arrival  time  of  the  signal  in  Chapter  3.  We  call  the  new  moments 
the  central  moments.  They  are  denoted  by  M^n)  and  are  given  by 


5- 


M(n)(z) 


(t-  W/0v )  <A*(z, t)A(z,t)>dt  -  l 


(1) 


M 


(0) 


n\  n--J 

M(j) 


j«0  m*'0^ 


n-0 ,1,2,  .... 


(58) 


We  note  ^  and  ^“0.  The  calculations  of  higher  order  central 

moments  follow  in  a  straightforward  manner  by  using  (58)  to  give 


M(2)(z)/M(0) 


T02  A0L  AA4 
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.  Temporal  Moments  and  Numerical  Results 


In  this  chapter,  we  give  the  physical  meanings  to  those  moments. 

The  first  moment  and  the  second  central  moment  have  the  physical  signifi¬ 
cance  as  the  mean  arrival  time  and  the  mean  square  width  of  the  received 
signal.  Then  the  skewness  and  kurtosls  are  defined  in  terms  of  the  second, 
third  and  fourth  central  moments.  These  quantities  are  used  to  describe 
the  rough  shape  of  the  signal.  In  order  to  get  more  concrete  understanding 
ve  assume  that  the  turbulence  has  the  power  spectrum  Introduced  by 
Shkarofsky.  Then  two  sets  of  parameters  are  used  to  obtain  the  numeri¬ 
cal  results.  These  parameters  describe  the  geometry  and  content  of  the 
medium.  All  of  them  are  stated  in  Section  3.1. 


3.1  Mean  Arrival  Time  and  Mean  Pulse  Width 


The  first  moment  and  second  central  moment  have  been  evaluated  and 
discussed  by  many  papers  in  the  past  [Mark,  1972;  Yeh  and  Yang,  1977c; 
Liu  and  Yeh,  1977;  Yeh  and  Liu,  1977b;  Yeh  and  Liu,  1979].  The  first 
moment  represents  the  time  position  of  the  energy  weighted  by  the  inten¬ 
sity  distribution.  Therefore  it  is  Just  the  arrival  time  of  the  signal. 
Here  we  rewrite  the  first  moment  expression  in  equation  (54)  and  denote 
t  as  the  arrival  time  of  the  signal.  So  we  have 

A? 

t  =  M(1)(z)/M(0)  =  z/v  -  -r1  (2Lz— Lz)X2Z-5/:  (6 

a  g  4c 

where  we  have  replaced  cZ1^2  by  v  ,  the  group  velocity.  Obviously  the 
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first  term  in  the  equation  above  represents  the  arrival  time  if  there 


were  no  random  irregularities  in  the  medium.  The  second  term  is  the 


excess  time  due  to  dispersion  and  random  scattering.  Here  we  have  to 
note  that  A?  is  always  negative,  then  the  excess  time  is  always  positive 


Since  Che  second  term  depends  on  X2,  Che  excess  time  decreases  rapidly 
as  Che  carrier  frequency  Increases. 


Nov  for  furcher  calculaCion,  we  lncroduce  a  power  spectrum  for  Che 
random  turbulence  which  was  presented  first  by  Shkarofsky  [1968] : 

<Koi0)(p‘3)/2*Q3  K  /,a0’£'z+<oT) 

$  ^  )  m  ■■■■'  ■  i  ■  -  .  •  £  (  ■  ^  i  y  ^  "  •  g  ^ 

C  (27t)3/2K(p_3)/2(k-0£0)  (£0^)p/2  5 


where  K  is  che  Hankel  function  of  imaginary  argument  and  kq-1/Lo.  The 
quantities  Lo  and  Iq  are  the  outer  and  inner  scales,  respectively,  o^2 
is  the  variance  of  the  relative  electron  density  fluctuation.  The  two- 
dimensional  correlation  function  corresponding  to  the  three  dimensional 
spectrum  defined  by  (63)  is  given  by 

2tt(k0i/o7+I72')  ^P”2^2Kf 

A,  (o) - - L- 

5  <o(<oAo)(P  )f K(p-3)/2(k°1o) 


When  A^(p)  is  expanded  as  in  equation  (37),  the  coefficient  A n«0,2,4, 
...  etc.,  can  be  easily  obtained  from  the  relation  A  ■  ^ **  I  _ 


n  n.  9pn  |o-0, 


yielding 


Aq  ■  /2it£q/<q  a^2  ^^p_2) ^(p-3) 


A,  -  -/-nco^‘0  V  lt(p-4)/2<'<0‘»)/K(ii.3)/2<.oio) 


A.  -  /.k0!/25J(J3  a.3  K(p.<i/2<,'0«0>/K(p.3)/2 


*  ■■  *27£o5  K(p-8)/2^<oi°^  /,K(p-3)/2<'K°£'0^ 

(68) 

A8  “  ^<0  /9*2  S,q  a^2 

(69) 


For  numerical  results,  we  take  Lq*10  kilometers  and  £o*10  meters 
and  cr^2«0.1.  Furthermore  we  use  two  sets  of  parameters  which  are  listed 
as  follows: 


Model  1 

Model  2 

Plasma  frequency 

f 

P 

10  MHz 

50  MHz 

Distance 

z 

500  km 

1000  km 

Width  of  random 
irregularity  slab 

L 

200  km 

500  km 

where  Model  1.  corresponds  to  conditions  that  can  occur  naturally  in 
the  equatorial  ionosphere  and  Model  2  corresponds  to  conditions  of  an 
ionosphere  disturbed  by  nuclear  explosions.  In  our  calculations  the 
carrier  frequency  f£  will  be  varied.  The  width  of  the  impressed  pulse 
is  To  and  is  also  varied  such  that  ^Tq  *100,  whic is  consistent  with 
the  assumption  made  in  the  derivation  of  the  temporal  moments  in  Section 
2.4. 

In  Figure  2  we  present  the  numerical  results  for  the  arrival  time. 

The  excess  time,  t  -z/v  ,  is  plotted  as  a  function  of  the  carrier  fre- 
a  g’ 

quency  for  both  Model  1  and  Model  2.  Except  in  the  region  fc<300  MHz 

in  Model  2,  these  two  straight  lines  in  this  figure  show  that  the  excess 

■>4 

times  are  proportional  to  f£  in  both  models.  The  little  digression  at 
low  carrier  frequency  end  in  Model  2  comes  from  the  failure  of  the  approxi¬ 


mation  Z;l. 


Now  let's  look  at  the  second  central  moment  given  by  equation  (59). 


It  has  been  identified  earlier  as  the  mean  square  pulse  width  of  the 
signal.  Denote  the  square  root  of  it  by  t.  So  we  have 


/q\  /|i\  To  AqL  4Ai* 

t2  3  MU'(z)/M(U'  -  —  +  X2Z‘3+  (3Lz2-3L2z+L3)X2Z-4 

1  ^ 

+  l-2'2-  T  >-3z  +  i  L-JX-Z-5 


It  is  easy  to  recognize  that  the  first  term  of  the  expression  above 

is  the  original  square  width.  The  other  terms  contribute  to  the  broadening 

of  the  signal  due  to  random  scattering  and  dispersion.  But  in  equation 

(70)  the  distortion  terms  that  come  from  the  high  order  dispersion  are 

absent.  The  reason  for  this  is  that  we  have  neglected  g2^2  term  and  all 

higher  ones  in  equation  (49).  Actually,  the  general  expressions  (47)  and 

(48)  contains  the  effects  due  to  higher  order  dispersion.  If  in  those 

terms  we  take  U22  term  in  equation  (49)  into  account  for  the  calculations 

of  and  an  extra  term  (z2/c2w  2To2)X2Z”3  must  be  added  to  the 

c 

right  hand  side  of  equation  (70).  This  quantity  describes  the  signal 
broadening  if  there  are  no  random  irregularities  in  the  medium.  Using 
the  numerical  values  in  either  Model  1  or  Model  2,  we  find  this  quantity 
is  much  smaller  than  any  distortion  term  in  equation  (70)  due  to  scattering 
and  hence  the  higher  order  dispersion  term  can  be  ignored. 

The  dependence  of  the  normalized  mean  pulse  width  t/(Tq/2)  on  the 
carrier  frequency  is  shown  in  Figure  3.  In  the  frequency  range  we  con¬ 
sider,  when  the  carrier  frequency  is  larger  than  500  MHz  for  Model  1 
and  that  is  larger  than  10  GHz  for  Model  2,  the  first  term  in  equation 
(70)  dominates  and  we  can  neglect  all  the  broadening  effects. 


In  Model  1,  on  the  other  hand,  if  the  carrier  frequency  is  between 
100  MHz  and  140  MHz,  the  term  containing  At,,  which  has  f  "***  dependence, 
is  most  important.  In  some  earlier  papers  [Sreenivasiah,  1976;  Leader, 
1979;  etc.]  A^(p)  in  equation  (37)  is  approximated  by  Aq+A2p2.  With 
this  approximation,  they  would  miss  the  A4  term  which  happens  to  be 
dominating  under  the  condition  above.  Then  their  results  will  be  in 
error  when  applied  to  our  problem.  In  Model  2,  when  the  carrier  fre¬ 
quency  fc  is  such  that  100  MHz  £  f£  _<  1  GHz,  the  last  term  in  equation 
(70)  is  most  important  and  we  can  approximate  t2  as 

A22 

t2  “  p-  (j  L2Z2-  -3“  L3Z+  i  L4)X4Z“5  (71) 
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which  has  f  dependence.  Note  that  under  this  condition,  the  normalized 
mean  pulse  width  t/(T0/2),  Is  proportional  to  f  since  in  ill  computa¬ 
tions  we  assume  u  Tn-100. 

c 

3.2  The  Third  Moment  and  Skewness 

Let  us  rewrite  the  third  central  moment  given  by  equation  (60)  as 

A2 

M(3)(z)/M(0)  -  XZ"5/2-  (2Lz-L2)X2(24+15X)Z-9/2 

c  c 


(12Lz3-18L2z2+12L3z-3L4)X2Z-1/2 


(12L2z3-22L3z2+15L4z-  —  L5)X4Z"13/2 


ir  (|  I-3-3-  \  ^  ‘■6>x6z-‘5/2 


The  first  term  on  the  right  hand  side  of  (7?)  comes  from  dispersion 
and  all  the  remaining  terms  from  random  scattering.  Note  that  A2  and  A., 


are  always  negative  and  A4  is  positive.  This  implies  that  the  third 
central  moment  is  always  positive.  Physically,  this  means  that  due 
to  propagation  effects  the  pulse  becomes  asymmetrical  with  a  trailing 
edge  longer  than  the  leading  edge.  If  under  certain  conditions  the 
third  central  moment  becomes  large,  the  received  pulse  must  be  then 
stretched  to  give  a  long  tall.  Therefore,  the  value  of  the  third 
central  moment  is  a  measure  of  asymmetry  of  the  pulse  relative  to  the 
arrival  time. 

In  Figure  A,  we  present  the  dependence  of  the  third  central  moment 
on  carrier  frequency  for  two  models.  In  Model  1,  when  the  carrier  fre¬ 
quency  f  is  larger  than  3  GHz,  the  first  term  in  equation  (72)  dominates, 
so  that 

M(3)(z)/M(0)  =  ^7  X  z-5/2 
c 

indicating  that  the  dispersion  is  more  important  and  the  third  moment 

~4 

has  £c  dependence.  When  f^l  GHz,  the  third  term  dominates, 

M(3)(z)/M(0)  =  •  -r  (12Lz3-18L2z2+12L3z-3L4)X2Z~u  ,2 
c 

indicating  that  the  random  scattering  is  more  important  and  this  moment 

—8 

is  proportional  to  f  .  In  Model  2,  when  the  carrier  frequency  is 

above  10  GHz,  the  term  owing  to  dispersion  effect  dominates  again. 

But  when  f  is  below  10  GHz,  no  term  can  always  dominate, 
c 

Let  Os  now  define  the  skewness,  s,  through  the  following  relation: 

M(3)(z)/M(0)  -  st3  (73) 

The  value  of  s  is  a  measure  of  the  extent  of  the  signal  asymmetry.  If  the 
signal  intensity  distribution  is  symmetric  about  the  arrival  time,  s  is 
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zero.  Roughly  speaking,  larger  s  means  bigger  relative  difference  of 
the  stretched  lengths  between  the  trailing  edge  and  the  leading  edge. 

In  Figure  5,  we  show  a  distribution  with  an  exponential  decay  in  trailing 
edge  and  a  Gaussian  decay  :.n  leading  edge.  This  distribution  has  skew¬ 
ness  s-1.03.  If  any  distribution  has  s  value  much  smaller  than  unity, 
we  may  conclude  that  any  of  the  following  three  possibilities  may  happen: 
1)  the  trailing  edge  decays  faster  than  the  exponential  decay;  2)  the 
leading  edge  decays  slower  than  the  Gaussian  decay;  3)  both  of  above. 

With  the  distribution  shown  in  Figure  5,  this  comparison  of  the  s  values 
will  give  us  a  coarse  idea  about  the  shape  of  the  received  signal. 

The  dependence  of  skewness  on  the  carrier  frequency  is  shown  in 
Figure  6  for  both  models.  In  Model  2,  when  the  carrier  frequency  f 
is  below  about  800  MHz,  s  has  nearly  a  constant  value  of  2.  When  f 

c 
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is  between  1.5  GHz  and  10  GHz,  s  is  proportional  to  f  and  when  f 

c  c 

is  above  15  GHz,  s  is  a  linear  function  of  f  .  In  Model  1,  except 

the  horizontal  part,  the  curve  is  roughly  parallel  to  that  of  Model  2. 

When  ffi  is  such  that  200  MHz^f  _j<2  GHz,  s  is  proportional  to  f and 

when  f  is  above  4  GHz,  s  is  inversely  proportional  to  f  again.  Note 
c  c 

that  if  the  carrier  frequency  is  above  500  MHz  for  Model  1  and  is  above 
3  GHz  for  Model  2,  the  values  of  s  are  much  smaller  than  1.  So  we  can 
say  in  those  frequency  ranges  the  received  signals  must  be  less  asymmetric 
than  the  distribution  shown  in  Figure  5,  if  by  some  other  means  we  can 
make  sure  that  the  signal  is  in  i  single  clump. 

3.3  The  Fourth  Moment  and  Kurtosls 

The  formula  of  the  fourth  central  moment  is  given  in  equation  (61). 
Note  that  since  every  term  in  the  right  hand  side  is  nonne.gati.ve,  the 
fourth  central  moment  is  nonnegatlve  as  it  should  be.  Again,  the  first 
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Figure  5.  A  distribution  with  an  exponential  decay  in  the 

trailing  edge  and  an  Gaussian  decay  in  the  leading 
edge.  It  has  skewness  1.03  and  kurtosis  2.94.  Note 
that  we  have  chosen  the  origin  of  the  horizontal  axis 
to  make  M  =0. 


Figure  6.  Skewness  s  as  6  function  of  the  carrier  frequency 
f  .  For  a  definition  of  s,  see  equation  (73)  in 
tfte  text. 
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term  is  the  original  fourth  central  moment  and  all  the  others  come  from 

the  propagation  effects  due  to  dispersion  and  random  scattering.  We 

plot  M^(z)/M^  as  a  function  of  carrier  frequency  in  Figure  7.  Quite 

similar  to  the  situation  of  the  second  central  moment,  when  f  is  above 

c 

1  GHz  in  Model  1  and  is  above  10  GHz  in  model  2,  the  first  term  in  equa¬ 
tion  (61)  dominates  and  we  can  omit  all  the  propagation  effects.  On 

the  other  hand,  in  Model  1  when  f  Is  below  500  MHz,  the  fourth  central 

c 

moment  is  proportional  to  f  and  it  can  be  approximated  by  the  single 
term  containing  Ag.  In  model  2,  when  the  carrier  frequency  is  between 
200  MHz  and  1  GHz,  this  moment  is  roughly  proportional  to  £  i.e., 

the  term  including  A24  is  moat  important. 

Now  we  define  the  kurtosis,  denoted  by  K,  from  the  following 
relation: 

M(4)(z)/M(0)  -  t4(K+3)  (74) 

where  x  is  the  mean  pulse  width  defined  before.  The  kurtosis  is  a  di¬ 
mensionless  measure  of  the  distribution  concentration  extent.  For  the 
distribution  shown  in  Figure  5,  the  value  of  K  is  2.94.  If  any  other 
one  clump  distribution  has  kurtosis  value  less  than  2.94,  it  must  be 
more  concentrated  than  the  distribution  shown  in  Figure  5. 

In  Figure  8,  we  show  the  dependences  of  K  on  the  carrier  frequency 
for  both  models.  Note  that  both  curves  have  the  same  shapes  as  the  cor¬ 
responding  ones  in  Figure  6.  Notice  also  that  K  is  always  larger  than  0, 
which  is  thr-  case  without  any  propagation  effect.  In  Model  2,  when  the 
carrier  frequency  is  below  3  GHz,  K  has  nearly  a  constant  value  of  30. 
When  the  carrier  frequency  la  above  700  MHz  for  Model  1  and  is  above 
4  GHz  for  Model  2,  the  values  of  kurtosis  are  less  than  2.94.  This 
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Figure  8.  Kurtosis  K  as  a  function  of  the  carrier  frequency 
f  .  For  a  definition  K,  see  equation  (74)  in  the 


implies  that  in  those  frequency  ranges  the  received  signal  intensity 
oust  be  more  concentrated.  So  we  can  conclude  that  in  those  frequency 
ranges  the  signal  must  decay  faster  than  the  exponential  decay  on  both 
sides  if  we  assume  the  received  signal  is  just  in  a  single  clump. 
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4.  Description  of  the  Signal  Tail 


In  many  applications  we  are  not  concerned  with  the  detailed  and 
exact  shape  of  a  propagating  signal,  such  as  in  binary  communication. 

The  important  thing  we  like  to  know  is  whether  an  error  will  be  committed 
in  the  decoding  process  owing  to  the  distortion  of  the  signal.  From  the 
results  In  previous  chapters,  we  know  the  distortions  due  to  propagation 
effects  can  be  described  in  two  separate  ways.  One  of  them  is  the  broadening 
of  a  pulse.  The  other  is  the  asymmetry  of  this  pulse,  which  when  trans¬ 
mitted  was  originally  symmetric.  In  this  latter  aspect,  we  find  the 
distorted  signal  always  has  a  longer  trailing  edge.  When  serious,  the 
received  signal  cannot  be  contained  in  its  original  bit,  it  will  extend 
its  energy  into  the  next  bit  or  even  the  third  bit.  This  may  produce 
errors  in  the  decision  making  process  for  the  neighboring  bits.  Since 
in  a  binary  communication  link,  we  mainly  worry  about  the  possibility 
of  committing  the  errors,  it  is  important  to  estimate  the  energy  quantity 
which  is  extended  outside  of  the  original  bit  itself,  especially  in  the 
trailing  edge. 

We  introduce  two  theorems  about  moments  first  in  this  chapter.  Then 
these  two  theorems  will  be  used  tc  give  an  upper  bound  for  the  signal 
energy  outside  some  time  interval  from  the  arrival  time  in  the  trailing 
edge.  The  parameter  values  of  those  two  models  in  the  last  chapter  will 
be  used  again  for  the  numerical  results. 

4.1  Two  Theorems  About  Moments 

Since  we  are  dealing  with  the  normalized  signal  intensity  distribu¬ 
tion,  this  distribution  density  as  a  function  of  time  is  always  positive 
and  the  total  distribution  can  be  normalized  to  unity.  These  properties 


satisfy  the  conditions  for  a  function  to  be  a  probability  distribution 
density  function.  So,  any  theorem  about  the  moments  of  the  probability 
distribution  as  a  function  of  random  variable  can  be  cited  in  our  re¬ 
search,  if  these  moments  are  defined  in  the  same  ways  as  the  temporal 
moments.  Here  let  us  Introduce  two  theorems  about  the  moments  of  the 
probability  distribution.  These  two  theorems  can  be  found  in  books  on 
statistics,  for  example,  the  book  written  by  Mises  [1964].  In  the  fol¬ 
lowing  statements,  we  will  not  distinguish  the  probability  distribution 
from  the  signal  Intensity  distribution,  since  the  results  have  no  dif¬ 
ference. 

Theorem  1:  If  two  distributions  with  distribution  density  function  I(x) 
and  I ' (x)  have  the  same  moments  up  to  the  nth  order,  the 
graphs  of  the  corresponding  cumulative  distribution  functions 
e(x)  and  e'(x)  «ust  have  at  least  n  "points  of  intersection". 
Note  that  the  cumulative  distribution  functions  and  the  distribution 
density  functions  are  related  by 

X 

e(x)  ■  |  I(x')  dx' 

Here  we  have  to  give  a  definition  to  "a  point  of  intersection".  We 
define  "a  point  of  intersection"  as  that  inside  some  closed  interval 
of  ;<  where  £(x)*e'(x)  and  just  outside  this  interval  e(x)-e'(x)  has 
opposite  signs  in  the  positive  side  and  negative  side. 

In  this  report,  we  will  not  try  to  prove  this  theorem.  But  from 
Theorem  1,  it  is  easy  to  derive  the  following  Theorem  2. 

Theorem  2:  Suppose  there  exists  a  distribution  l(x)  which  has  an 
ra-step  increasing  cumulative  distribution  function  e(x) 
and  has  first  2m  order  moments  M^,  M^, 


*  *  M 


— ( 2m-l ) 


If  I ' (x)  is  any  other  distribution  with  the  same  2m  moments 
as  above,  then  the  graph  of  its  cumulative  distribution 
function  e'(x)  passes  through  each  "step"  and  each  "riser" 
of  this  o-step  function  e(x)» 

As  an  example,  ve  show  a  3-step  case  in  Figure  9.  In  this  figure,  for 
the  purpose  of  our  application  we  have  replaced  the  variable  x  by  t  and 
denote  the  positions  of  the  "risers"  by  ti,  t2,  t3  and  the  heights  rr 
the  respective  riser  by  e j ,  £2,  £3.  Note  that  £1,  £2,  £3  must  be  posi¬ 
tive  and  the  moments  up  to  fifth  order  are  concerned.  Since  in  the 
statements  of  Theorem  2,  we  have  used  the  central  moments  to  specify 
the  m-step  increasing  function,  the  origin  of  the  horizontal  axis  is 
the  arrival  time  of  the  signal. 

But  for  given  moments  ,  ....,  M^2™  such  an 

m-step  increasing  cumulative  distribution  function  can  exist  only 
when  some  conditions  are  satisfied.  For  the  case  m-3,  these  conditions 
are 


M 


(0) 


>  0, 


M(0)  M(1> 

m(1)m{2) 


M<°>  M(1>  M(2) 
M(1>  M(2)  M(3) 

m(2)  m(3)  m<4) 


>  0 


4.2  Application  of  Those  Theorems 

In  the  case  of  a  3-step  function  stated  above,  for  the  purpose  of 
calculating  the  positions  and  heights  of  "risers",  i.e.,  t \ ,  t2,  t3, 
el»  c2»  £3  (b  unknowns)  in  Figure  9,  the  quantities  of  moments  up  to 
fifth  order  (6  knowns)  are  needed.  But,  if  we  have  the  information  ^ 
only  up  to  fourth  moment,  we  still  can  construct  a  3-step  increasing 
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Figure  9.  A  3-step  function  and  a  increasing  curve  whose 

corresponding  distribution  density  functions  have 
the  same  moments  up  to  fifth  order,  c.  represents 
the  height  of  the  "riser"  at  t^  i«lt2,3.  Note 
that  ei+e2+£3"l.  In  this  case,  these  two  curves 
have  at  least  five  "points  of  intersection".  Since 
the  height  between  point  P  and  Q  is  the  fractional 
signal  energy  beyond  t3,  £3  is  an  upper  bound  of  this 
fractional  energy. 
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function  except  that  any  one  of  tj ,  t2*  t3,  ei,  t2»  £3,  must  be  left 
undetermined.  This  undetermined  one  can  be  chosen  arbitrarily  within 
some  restricted  ranges. 

The  conditions  under  which  such  a  3-step  increasing  function  can 
exist  have  been  listed  at  the  end  of  Section  4.1.  They  can  be  satisfied 
if  the  cumulative  distribution  function  of  our  signal  intensity  distri¬ 
bution  id  increasing  (not  constant)  at  least  at  3  points.  This  is  true, 
since  the  signal  Intensity  is  monotonlcally  increasing  at  least  in  some 
range  of  time.  The  proof  for  the  satisfaction  of  those  conditions  can 
be  seen  just  by  reviewing  the  process  of  deriving  those  conditions, 
which  we  omitted. 

Since  we  are  trying  to  estimate  the  signal  energy  outside  some 
time  interval  from  the  arrival  time  in  the  trailing  edge,  among  t^,  t2» 

1 3 ,  ei,  e 2 •  e 3  we  will  let  t3  unfixed.  Now  let  us  compute  tj ,  t2,  cj,  t2» 
e3  in  terms  of  t3,  T,  s,  K.  The  starting  point  of  the  calculation  is 
the  equivalence  of  those  moments  between  our  signal  intensity  distribu¬ 
tion  and  the  distribution  with  3-step  increasing  cumulative  distribution 
function  e(t).  The  corresponding  distribution  density  function  I ( t )  of 
e(t)  can  be  expressed  by  a  sum  of  delta  functions  as 

I(t)  -  Ei6(t-ti)+e26(t-t2)+e 3o(t-t3)  (75) 


Then  we  have  relations  as  follows: 


M(o)/M(°)  „  ei+e:+E? 


M(l)/M(u)  .  £  it:+t:t:+<.  it 


(76) 


(77) 
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M(2)/M(0)  -  eit12+e2t22+£3t32 


M^/M^  •  c 1 1 1 ^+e2t2  ^+c 3 1 3  2 


M(4)/M<0)  -  eit14+e2t24+e3t34 
-  t4(K+3) 


From  equation  (75),  (76)  and  (77),  it  is  easy  to  solve  ej,  e2,  e3  in 


terms  of  ti,  t2,  t3.  Th«  results  are 


xz+t2t3 

( t 3— t i ) (t2-ti) 


T  1  ^  t  3 

(t3~t2) (tl“t2) 


T4’+tlt2 

e  3  ■  -  (8 

(t!-t3)(t2-t3) 

Substituting  these  results  in  equations  (79)  and  (80),  we  can  obtain  the 
expressions  for  t i ,  t2  in  terms  of  t3,  x,  s  and  K.  They  are  as  follows: 


'•Fr 


where 


T2t3(KH-2)4-T3S-t32ST 

T2+t3TS-t32 


T2  [T‘’-S2-t3TS-T2(K+3)+t32] 


q  -  - - - 

T2+t  3  TS  — 1 3  2 

Now  we  have  constructed  the  3-step  Increasing  function  e(t).  Since 
we  know  the  cumulative  distribution  function  e'(t)  of  our  signal  intensity 
distribution  is  always  Increasing  and  must  huve  at  least  5  "points  of 
intersection"  with  e(t),  the  graph  of  e'(t)  must  pass  through  the  third 
"riser"  of  e(t)  at  t3  (see  Figure  9).  Furthermore,  since  the  height 
from  this  Intersection  point  to  the  top  step  is  the  fractional  signal 
energy  beyond  t3,  £3  is  just  an  upper  bound  of  this  fractional  en  ,y. 

In  a  binary  communication,  we  can  always  put  t3  at  the  boundary  of  a  bit 
and  then  compute  the  fractional  signal  energy  extended  into  the  neighboring 
bits  in  the  trailing  edge.  But,  we  have  to  remember,  £3  can  be  chosen 
only  within  some  ranges.  These  ranges  can  be  found  from  the  conditions 
which  require  that  ti  and  t2  must  be  real  and  e^,  eo,  c3  must  be  real 
and  positive. 

4.3  Numerical  Results 

In  this  section,  we  use  the  model  parameters  given  in  Section  3.1 
once  again.  We  plot  the  curves  of  £3  as  a  function  of  the  carrier  fre¬ 
quency  in  Figure  10.  Two  values  of  t3  are  assigned,  one  at  Tg,  the  other 
et  2Tq .  The  shapes  of  those  four  curves  are  quite  similar.  In  both  models, 
£3  asymptotes  to  about  0.1  for  t3*Tg  and  to  about  0.005  for  t3*2T0.  Also 
we  find  when  the  carrier  frequency  is  larger  than  1  GHz  In  Model  1  and 
larger  than  10  GHz  in  Model  2,  £3  is  approximately  equal  to  the  individual 
asymptotic  value.  These  two  values  of  the  carrier  frequency  coincide 
roughly  with  those  beyond  which  the  propagation  effects  on  tin-  pul'.*- 


width  can  be  neglected,  as  shown  in  Figure  3.  Also  from  Figure  6, 
when  the  carrier  frequencies  are  above  those  two  values  in  two  models, 
respectively,  skewness  values  are  less  than  0.001,  which  is  much  smaller 
than  that  of  the  distribution  in  Figure  5.  So  we  can  conclude  that  the 
signal  shape  is  nearly  unchanged  for  the  carrier  frequency  range  above 
those  two  values  in  two  models,  respectively,  Really,  for  the  undis¬ 
torted  signal,  the  exact  fractional  energy  beyond  Tq  is  0.023  and  tht.t 
beyond  2Tq  is  5xl0~^,  which  naturally  are  smaller  than  their  upper  bound 
0.1  and  0.005  given  above,  respectively.  Because  we  have  the  restrictions 
on  the  choice  of  t3,  we  can  not  complete  those  curves  in  Figure  10  when 
the  carrier  frequency  is  below  some  value. 


fc  (GHz) 

Figure  10.  An  upper  bou^.c  of  the  fractional  energy  beyond 
1 3 »  Ei  *  function  of  the  carrier  frequency. 
Two  values  of  1 3  are  assigned. 


Summary  and  Conclusion 


Starting  from  the  Helmholtz  wave  equation  (10),  we  derived  an 
equation  (17)  for  the  two-frequency  two-position  mutual  coherence 
function  F(p).  In  order  to  compute  I0,  r 1 •  T2,  F3,  r4,  a  set  of 
equations  for  and  W^'  were  derived.  r0,  F 1 ,  T2,  F3,  r4  were  eval¬ 
uated  and  then  temporal  momenta  and  temporal  central  moments  up  to  the 
fourth  order  were  computed.  Two  models  for  the  geometry  and  ionospheric 
parameters  were  used  to  obtain  numerical  results. 

We  next  consider  a  narrow-band  Gaussian  envelope  carrier  signal 
being  impressed  at  z™0.  After  propagating  through  a  turbulent  plasma, 
owing  to  dispersion  and  random  scattering,  this  originally  symmetric 
signal  is  broadened  and  becomes  asymmetric.  The  trailing  edge  is  longer 
than  the  leading  edge.  From  the  information  we  obtained,  we  can  not  tell 
whether  the  received  signal  is  just  in  a  single  clump.  But  if  we  can  make 
sure  of  it  by  some  other  means,  we  compare  the  skewness  and  kurtosis  values 
of  the  average  signal  intensity  distribution  with  those  of  the  distribution 
shown  in  Figure  5  and  then  get  a  rough  iden  about  the  shape  of  the  re¬ 
ceived  signal. 

The  results  also  showed  the  dominating  propagation  effect  between 
dispersion  and  random  scattering.  As  the  arrival  time  and  pulse  width 
are  concerned,  the  random  scat.  ;•  ing  effect  is  more  important,  since  the 
propagation  effect  will  * e  much  smaller  if  there  are  no  random  irregu¬ 
larities  in  the  medium.  As  far  as  the  signal  asymmetry  is  concerned, 
the  propagation  effect  mainly  comes  from  dispersion  in  high  frequency 
part  and  from  random  scattering  in  low  frequency  part  in  the  carrier 
frequency  range  we  considered.  Finally,  for  the  extent  of  signal  con¬ 
centration,  random  scattering  effect  is  dominating. 


Although  we  can  not  know  the  exact  shape  of  the  received  signal 
intensity  distribution,  we  can  find  an  upper  bound  for  the  fractional 
signal  energy  beyond  some  time  distance  from  the  arrival  time.  In  a 
binary  communication,  this  information  may  help  us  to  predict  the  errors 
in  a  decoding  process.  But  to  do  it,  we  need  more  investigations. 


References 


Leader,  J.  C.  (1978),  Intensity  fluctuations  resulting  from  partially 
coherent  light  propagating  through  atmospheric  turbulence,  J.  Opt. 
Soc,  Am.,  69,  pp. 73-84. 

Liu,  C.  H.,  and  K.  C.  Yeh  (1977),  Propagation  of  pulsed  beam  waves 
through  turbulence,  cloud,  rain,  or  fog,  J.  Opt.  Soc.  Am.,  67, 

No.  9,  September. 

Liu,  C.  H.,  and  K.  C.  Yeh  (1978),  Pulse  propagation  In  random  media, 

IEEE  Trans,  on  Antennas  and  Propagation,  Vol.  AP-26,  No.  4,  July. 

Mark,  W.  D.  (1972),  Characterization  of  stochastic  transients  and  trans¬ 
mission  media:  The  method  of  power-moments  spectra,  Journal  of 
Sound  and  Vibration,  22(3),  pp. 249-295. 

Mises,  R.  v.  (196^),  Mathematical  Theory  of  Probability  and  Statistics, 
Sec.  4  and  5,  Chapter  VIII,  pp. 384-396,  Academic  Press. 

Shkarofsky,  I.  P.  (1968),  Generalized  turbulence  space-correlation  and 

wavu-number  spectrum- function  pairs,  Can.  J.  Phys.,  Vol.  46,  pp.2133- 
2153. 

Sreenivasia'h,  l.,  A.  Ishimaru,  and  S.  T.  Hong  (1976),  Two-frequency  mutual 
coherence  function  and  pulse  propagation  in  a  random  medium:  An 
analytic  solution  to  the  plane  wave  case.  Radio  Sci.,  11,  775-778. 

Tatarskii,  V.  I.  (1.971),  The  Effects  of  the  Turbulent  Atmosphere  on  Wave 
Propagation,  U.S.  Department  of  Commerce,  National  Technical  Infor¬ 
mation  Service,  Springfield,  Va , 

Yeh,  K.  C.t  and  C.  H.  Liu  (1977a),  An  investigation  of  temporal  moments 
of  stochastic  waves,  Radio  Sci.,  Vol.  12,  No.  5,  pp. 671-680,  Sept.- 
Oct. 


&A 


Yeh,  K.  C.,  and  C.  H.  Liu  (1977b),  Pulse  delay  and  pulse  distortion  by 

random  scattering  in  the  ionosphere,  AGARD-CP-244,  Aspects  of  Electro¬ 
magnetic  Wave  Scattering  in  Radio  Communication,  edited  by  A.  N. 

Ince. 

Yeh,  K.  C.,  and  C.  C.  Yang  (1977c),  Mean  arrival  time  and  mean  pulsewidth 
of  signals  propagation  through  a  dispersive  and  random  medium, 

IEEE  Trans,  on  Antennas  and  Propagation,  Vol.  AP-25,  No.  5,  Sept. 

Yeh,  K.  C.,  and  C.  H.  Liu  (1979),  Ionospheric  effects  on  radio  communica¬ 
tion  and  ranging  pulses,  IEEE  Trans,  on  Antennas  and  Propagation, 

Vol.  AP-27 ,  No.  6,  Nov. 


Appendix  B 


A  DETERMINISTIC  STUDY  OF  PULSE  PROPAGATION  IN  AN 
ELECTRON  BUBBLE  MEDIUM 


BY 

MICHAEL  ROBERT  TUCKER 
B.S.,  University  of  Illinois,  1980 


THESIS 

Submitted  in  partial  fulfillment  of  the  requirements 
for  the  degree  of  Master  of  Science  in  Electrical  Engineering 
in  the  Graduate  College  of  the 
University  of  Illinois  at  Urbana-Champaign,  1981 


Urbana,  Illinois 


97 


UNIVERSITY  OF  ILLINOIS  AT  URBANA-CHAMPAIGN 


THE  GRADUATE  COLLEGE 


July,  1981 


WE  HEREBY  RECOMMEND  THAT  THE  THESIS  BY 


_ MICHAEL  ROBERT,  TUCKES _ 

fvtttt  rn  A  DETERMINISTIC  STUDY  OF  PULSE  PROPAGATION  IN  AM _ 

_ ELECTRON  BUBBLE  MEDIUM _ 

BE  ACCEPTED  IN'  PARTIAL  FULFILLMENT  OF  THE  REQUIREMENTS  FOR 

THE  DEGREE  OF _ HASTES.  OF  SCIENCE _ 


U- 

/ 


Director  of  Theji*  Research 


Head 


Id  of  D/i 


partment 


Committee  on  Final  Examination! 


Chairman 


r  Required  tor  loctor'j  decree  but  not  tor  marter’i. 


RR 


ACKNOWLEDGMENT 


I  would  like  to  thank  Dr.  K.  C.  Yeh  for  his  guidance  throughout 
my  year  of  research.  I  would  also  like  to  express  my  gratitude  to  Mrs. 
Linda  Houston  for  typing  this  thesis. 

This  work  was  supported  by  the  National  Science  Foundation  through 
Grant  No.  ATM-80-07039  and  by  the  Air  Force  Geophysics  Laboratory 
through  contract  F19623-78-C-0195. 


Appendix  B 


Table  of  Contents 

Page 

1.  Introduction . 

2.  The  Parabolic  Equation .  92 

2.1  Derivation. . gp 

2.2  Determination  of  £j(r) . 95 

3 .  The  Electron  Bubble .  93 

3.1  The  Bubble  Model. .  g8 

3.2  Program  SHEET  1 .  100 

4.  Moments .  108 

4.1  Definition . . .  \  qq 

4.2  Physical  Significance  . 

5 .  Da  ta . . .  j  i  ] 

6.  Analysis . 

References . . . . . .  . .  144 

Appendix  A .  146 

Appendix  B. .  171 


90 


1. 


Introduction 


As  a  radio  wave  propagates  through  the  ionosphere,  random  varia¬ 
tions  in  the  electron  density  will  cause  the  wave  to  fluctuate.  These 
fluctuations  are  called  scintillations.  Scintillations  are  especially 
apparent  in  regions  called  electron  bubbles.  Electron  bubbles  are 
regions  in  the  ionosphere  where  electron  density  is  very  low.  Also 
very  sharp  electron  density  gradients  are  present. 

In  this  report  I  will  make  a  deterministic  study  of  pulse  propa¬ 
gation  through  a  model  electron  bubble  at  several  different  carrier 
frequencies.  Because  an  electron  bubble  is  a  dispersive  medium,  one 
would  expect  a  pulse  propagating  through  a  bubble  to  be  distorted. 

I  will  obtain  a  graphical  relationship  between  the  carrier  frequency 
and  the  pulse  distortion.  To  obtain  this  relationship  I  will  calcu¬ 
late  the  moments  of  the  pulse  envelope  after  the  pulse  has  been  trans¬ 
mitted  through  an  electron  bubble. 


1 


\ 


2.  The  Parabolic  Equation 

Microwave  propagation  through  an  electron  bubble  is  a  strong 
fluctuation  problem  because  the  log  amplitude  variance  is  greater 
than  (.2-. 5).  The  parabolic  equation  method  [Tatarskii,  1971]  will 
be  used  to  study  this  problem. 

Figure  1  shows  the  geometry  of  the  problem.  A  wave  propagates 
from  z»0  to  z»Zq  through  an  electron  bubble.  In  this  report,  the 
geometry  is  considered  to  be  two  dimensional  so  that  there  is  no  y 
dependence. 


2.1  Derivation 

The  dielectric  permittivity  at  any  point  r-(z,x)  is  written  as 
an  average  or  background  value  plus  a  fluctuating  part. 

er<r)  "  <ef  (*) > f  1-t-c i  (jr>  ] 

The  average  wave  number  is  defined  to  be 


k  2  ■  W  2  U  q  C  q  *:  e  ( z )  > 

k2  -  k02<cr(z)> 


(1) 


(2) 

(3) 


The  source  free  Maxwell's  Equations  are  the  starting  point  for  the 
derivation  of  the  parabolic  equation. 


V*E(r)  -  -iu>u  oH(r) 


7<H(r)  “  i»e  (r)E(r) 


(4) 

(5) 


Taking  the  curl  of  (A)  gives 


7>7*E(r)  -  -lwij0(7*H(r)) 


(6) 


9? 


Substituting  (5)  into  (6)  results  in 


V*VxE.(r)  "  -iwyo(iw£r (r)E(r_))  (7) 

Using  (3)  and  the  vector  identity  7xV*A-grad  div  A-72A  gives 

grad  div  E(r)  -72E(_r)  -  k0 2er (r)E(r)  (8) 

If  the  irregularities  are  large  compared  to  the  signal  wavelength 
then  div  E(r)  -  0  [Tatarskii,  1971].  Equation  (8)  reduces  to 

V2E(r)+k02£r(r)E(r)  -  0  (9) 

The  x  component  of  the  electric  field  E(r.)  can  be  written  as 

Ex(r)  ■  U(r)eikz  (10) 

where  e  iu,t  dependence  is  understood,  propagation  is  along  the  z  direc¬ 
tion  and  U(_r)  is  complex.  From  (1) 

k02er(£)  *  k2 (1+e i (r) )  (11) 

Substituting  (10)  and  (11)  into  (9)  gives 

V2U(r)eikz+k2(l+£1(r))U(r)eikz  -  0  (12) 


or 


ikz  ,  3 


W-  +  £*  ♦  &  »<£> 


ikz 


3y 


3z‘ 


+k2U(ir)eikz+k2£i  (r)U(r)  -  0 


(13) 


also 


U(r)e 


ikz  „  ..ikz  ii 


ikz  3 


U(£)+2ike  ^  U(r)-U(r)k‘e 


ikz 


(14) 


9/1 


Substituting  (14)  into  (13)  gives 


32  >  ikz  8*  ...  .  ikz  ikz  3‘ 

Vl(£)'  * 


U(r)+2ike  U(r) 


-k2U(r)elk*+k2U(r)®lk*+k2e1(r)U(r)«lk*  -  0 


Simplifying  (15)  results  in 


V2U(r)+2ik  1“  U(r)+k2£l(r)U(r) 


The  scale  size  l  of  the  medium  is  defined  to  be  the  average  distance 
over  which  the  fluctuating  part  of  the  dielectric  permittivity  remains 
correlated.  If  l>>\  then 


[kfj  U(r)|  » 


[Tat;  rskii,  1971]  .  So  V2  can  be  replaced  by  the  transverse  Laplacian 


T  3x2  9y- 


and  (16)  becomes  the  parabolic  equation. 


VT2U(r)+2ik  ~  U(r)+k2Gl(r)U(r)  -  0 


The  strong  fluctuation  problem  is  reduced  to  solving  the  parabolic 
equation  for  a  medium  with  ei(_r)  variations.  These  variations  can 
be  treated  stochastically  where  ei(r)  is  a  random  process.  However 
in  an  electron  bubble  medium,  very  sharp  electron  density  gradients 
occur.  As  will  be  seen  in  the  next  section,  these  will  lead  to  large 
fluctuations  in  ei(r_).  In  this  case  e  i  (O  will  be  treated  determinis¬ 


tically. 


2.2  Determination  of  ei(r) 
Recalling  Equation  (1) 


er(r)  -  <er(z)>[l+e1(r)] 

where : 

er(£>  is  the  relative  dielectric  permittivity  at  any  point  r 
<s  (z)>  is  the  background  relative  dielectric  permittivity  for 
a  height  z  (the  value  at  the  edge  of  the  bubble). 
ei(_r)  is  the  varying  part  of  the  relative  permittivity. 

The  background  relative  dielectric  permittivity  is  given  in  terms  of 
the  background  plasma  frequency. 


<er(z)>  •  1  - 


w  2(z)  - 


N0(z)e2 


Nq(z)  is  the  background  electron  density  for  a  height  z. 
e  is  the  electron  charge, 
m  is  the  electron  mass. 

£0  is  the  free  space  dielectric  permittivity. 

The  relative  dielectric  permittivity  at  any  point  r_  is  given  by 

u  2(r) 


er(r)  -  1  - 


w  2  (r) 

P  - 


N(r)e2 


where  N(r)  is  the  electron  density  at  any  point  _r 


N(r)  -  N0 ( z) -AN (r) 
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(25) 


Substituting  (27)  and  (21)  into  (28)  and  solving  for  ei(r)  results  in 


ei  (r) 


AN(r) 

C  N0(z) 


]  /  (1- 


(29) 


So  the  fluctuating  part  of  the  relative  dielectric  permittivity  is 
dependent  on  the  frequency  of  the  wave,  the  background  plasma  fre¬ 
quency  and  the  fluctuations  in  the  electron  density.  In  an  electron 
bubble,  the  electron  density  fluctuations  will  be  spiky,  thus  |--x(r_) 
cannot  be  treated  as  a  statistically  homogeneous  random  process.  So  the 
statistics  ofs^(r)  cannot  be  uniquely  determined  by  its  mean  and  auto¬ 
correlation  (or  by  its  power  spectrum).  Thus  a  deterministic  model  is 
used  to  study  wave  propagation  in  a  medium  of  this  type. 


f)'7 


3. 


The  Electron  Bubble 


3.1  The  Bubble  Model 

In  several  observation  programs  near  the  magnetic  equator,  scintil¬ 
lation  at  Gigahertz  frequencies  has  been  related  to  the  development  of 
electron  bubbles.  Soon  after  sunset,  thin  layers  of  irregularities 
develop  below  the  F  layer  and  quickly  start  to  rise.  The  irregularity 
layer  thickens  and  develops  into  a  region  of  depleted  electron  density. 
These  regions  are  called  electron  bubbles. 

In-situ  rocket  and  satellite  measurements  are  used  to  provide 
further  information  about  electron  bubbles.  Specifically,  in-situ 
data  measured  by  McClure  et  al.  (1977)  was  used  to  develop  the  model 
of  the  bubble  used  in  this  report.  The  data  is  shown  in  Figure  2. 

Figure  2  shows  electron  density  measurements  along  a  satellite  path 
and  clearly  indicates  irregularities  and  sharp  electron  gradients. 

I  will  use  the  model  developed  by  Wernik,  (1979).  Wernik  in  modeling 
the  bubble  assumed  a  two  dimensional  model.  Also  he  a1  3umed  the 
horizontal  variations  of  the  electron  density  in  the  bubble  to  be  that 
given  in  Figure  2  at  every  height  z.  The  background  density  (the  den¬ 
sity  of  the  background  ionosphere)  is  assumed  to  follow  a  parabolic 
profile . 

N0(z)  -  Nm  (Wz-Z^/Ho^G  (30) 

where  is  the  maximum  electron  concentration  at  height  Z^,  and  Hq 
is  the  scale  height.  During  the  course  of  bubble  development,  its 
leading  edge  becomes  sharper  while  its  trailing  (lowe.)  edge  becomes 
more  blunt.  Wernik  defined  two  stages  of  development,  the  initial 
stage  and  the  developed  stage.  G  is  a  weighting  function  which  takes 


the  different  stages  into  account.  In  this  thesis  I  will  use  the 
developed  stage.  The  weighting  function  G  is  shown  in  Figure  3  for 
the  developed  stage. 

To  calculate  the  electron  density  at  a  point  r_»(z,x)  inside  the 
bubble,  the  following  equation  is  used 

AN(r) 

N(r)  -  N0(z)  -  [  ]  N0 (z)  (31) 

where  Nq(z)  is  given  by  (30)  and  AN(r)/No7z)  is  obtained  from  Figure  2. 

In  order  to  obtain  the  solution  to  the  parabolic  equation,  the 
fluctuating  part  of  the  relative  dielectric  permittivity  must  be  solved 
according  to  Equation  (29).  When  f,  the  frequency  of  the  wave  is  much 
larger  than  the  plasma  frequency  f^,  which  is  always  the  case  in  this 
thesis,  then  (29)  can  be  simplified  to 

8O.6(10~6)AN(r) 

£l(r)  «  - j? - ~  (32) 

where  AN(r)*N0(z)-N(r)  and  is  in  electrons/cm3  and  f  is  in  MHz. 

3.2  Program  SHEET  1 

Wernik  (1979)  developed  a  program  called  SHEET  to  model  wave 
propagation  through  ar  electron  bubble.  I  have  rewritten  and  revised 
his  program  to  mike  it  run  more  efficiently.  My  program  is  called 
SHEET  1,  SHEET  1  solves  the  wave  propagation  problem  In  the  electron 
bubble  model  described  in  the  previous  section  by  obtaining  solutions 
for  the  parabolic  equation.  The  following  analytic  solution  to  the 
parabolic  equation  was  used  by  Wernik  in  his  SHEET  program  and  is  used 
in  ray  SHEET  1  urogram. 


loo 


The  parabolic  equation  (19)  was  obtained  in  section  2.1. 


7T''U(£)-2ik  ~  U(r)  +  k2£l(r)U(r)  -  0 
Since  the  bubble  model  used  in  this  report  is  2  dimensional,  r*(z,x)  so 


VT2  -  S2/3x2  (33) 

thus  (19)  reduces  to 

-2ik  |j  U(r)  +  U(r)+k2ei(r)U(r)  -  0  (34) 

The  following  substitutions  are  made 

t  »  kz  b  ■  kx  (35) 

Substituting  (35)  into  (34)  gives 

.L.u(t,b)+-|  +  1  e1(t,b)U(t,b)  -  0  (36) 


The  entire  medium  la  partitioned  into  a  grid  as  shown  in  Figure  4. 
The  step  sizes  in  (z,x)  space  are  P  and  M.  These  must  be  changed 
into  step  sizes  in  (t,b)  space  by 


t  -  P  •  k 
H  -  M  •  k 

To  solve  (36)  numerically,  the  implicit  Crank-Nicholson  difference 
scheme  is  used.  This  scheme  uses  the  following  analog  to  the  second 
derivative 


(37) 


-^U(t.b)  -  jjy  lo(ll£j[  - 


2U 


,n+l 


°£{> 


(1kj)(uJ+1-2uJ  +  U"_1)]  (38) 


10? 


and  che  following  analog  to  the  first  derivative 


U(t,b)  -  (Uj+1  -  u")/t  (39) 

where 


j  ■  0,1,2, ... ,  J 
n  *  0,1,2,...,  N 
o  ■  1/2 


J  is  the  number  of  steps  in  the  b  direction  (also  in  the  x  direction) 
and  N  is  the  number  of  steps  in  the  t  direction  (also  in  the  z  direction). 
The  term  ej (t,b)U(t,b)  is  expressed  as 

n+1  n 

— 1 - [aUj+l  +  (1~0)Uj]  (40) 

The  f  function  is  defined  to  be 

f"  -  +®l  j)/2  (41) 

Substituting  (41),  (40),  (39)  and  (38)  into  (36)  gives 


±Vun+1  + 

2tF  j+1 


/ ,  it 

(1  -  r5  + 


it  o 
2 


fn)(Jn+l 


it 

2H" 


aU 


n+1 

j-1 


it 

tiT 


-(l-a)uj+1+[l+ 


-  ^(l-7)fjiuj  -  (l-o)  Uj_x 


(42) 


A  problem  occurs  at  the  boundaries  of  the  bubble  where  the  equation 
depend  on  Uq  and  "False"  boundries  are  set  up  with  UrJ-O  and 

Uj+^*0.  The  effect  of  these  bounories  does  not  penetrate  deeply  into 
the  field. 


J  0/1 


Equation  (42)  can  be  expressed  as  a  matrix  multiplication 


BnUn+1  -  (Bn-x-An)Un  n-0,1, .  . .  ,N 


(43) 


so  for  each  step  in  the  z  direction  there  is  a  matrix  A  and  a 
matrix  Bn. 


nfi 

Uo 

• 

• 

• 

e 

• 

• 

• 

■1 

Q11  A11 

B  -t*A 

uS 

. 

• 

• 

•5 

(44) 


The  matrices  An  and  Bn  are  tridiagonal.  Only  the  main  diagonal  and 

the  diagonal  on  each  side  of  the  main  diagonal  have  non-zero  elements. 

The  reason  for  this  is  because  (42)  depends  only  on  U?*1,  u"+?',  U1?,  u” 

J  J~l  J  j+l 

and  The  elements  of  the  matrices  ere  found  from  (42) 


_n  ito  .  .  £n 

Bj,j  "  1  -H*-+ 


(45) 


Rn  m  -n  ixo 

J+l.J  J.J+1  2HF 


(46) 


.n  -i  .  i  _n 

AJ  »j  "  H2  +  2  fj 


(47) 


.  n  .n 

J+l,j  “  J.J+l  “  H‘ 


(48) 


for  j-0,1,2 . . J 


,n+l 


Starting  from  (43)  U  is  solved 
(Bn-x*An)Ur 


bV4*1  -  _ * 


U1*1  -  (l-(Bn)"1TAn]Un 


(49) 


I  Ob 


un+1  .  Un  -  TXn 


(50) 


where 


Xn  -  (Bn)_1  AnUn 


(51) 


or 


bV*  -  aV1 


(52) 


Equation  (52)  is  solved  in  SHEET  1  for  X  using  Gaussian  elimination 
via  the  Thomas  algorithm  [Ames  (1969)].  Since  Bn  is  tridiagonal, 
Equation  (52)  can  be  written  as 


cciXi4-bbiX2  ■  dj 


(53) 


aaJXJ~l+CCJXJ  *  dJ 


where: 


aSj ,  ccj  and  bb^  are  elements  of  matrix  B  ,  and 


x0  -  xJ+1 


dJ  *  Vl-1  +  +  bj°J+l 


di  -  citj"  +  biUo 


(54) 


dJ  *  ajUS-l  +  C7UJ 


,n_„n 


where : 

a j ,  Cj  and  b^  are  elements  of  matrix  \n ,  and  U'o*U“+^-0. 

The  Thomas  Algorithm  is  used  to  solve  (52).  A  solution  is 
assumed 

*j  '  Vl*.)+l+£>! 

lor. 


(55) 


where: 


since 


vi  ■ 


Vi  ■  ^j~aaj®^^aajaj+ccj) 


a2  ”  -bbi/cci  $2  *  d2/cci 


"j+i  '  0 


*J  •  so+l 


*j  ■  “fjV""/'/"]' 


After  the  Xn  matrix  is  solved,  (50)  is  used  to  obtain  the  solution 


of  the  parabolic  equation  at  the  next  t  (or  z)  step. 


Moments 


u 

4.1  Definition 

In  this  thesia,  pulses  of  three  different  carrier  frequencies 
and  two  different  pulse  widths  are  propagated  through  an  electron 
bubble  mediuto.  In  order  to  analyze  the  output  pulse,  the  method  of 
temporal  moments  will  be  used.  Specifically,  the  first  four  temporal 
moments  will  be  used  to  giv  ■*  a  description  of  the  output  pulse  envelope. 
The  pulse  incident  on  the  electron  bubble  can  be  described  as 

ao 

Pi(£,u)  *  ' *d(jj  (62) 

— OO 

Where  F(uj)  Is  the  frequency  spectrum  of  the  pulse.  Assuming  propagation 
parallel  to  the  z  axis,  the  pulse  at  a  receiver  at  (zq,Xo)  is 

QO 

p0(zo,x3,t)  *  F(u))U(zo»Xf\,w)e~*  ^t-kz^dw  (63) 

—CO 

where  U(zq,x0,u)  is  the  parabolic  equation  solution.  The  output  pulse 
can  alternately  be  described  as 

io)  (t-z/c) 

Po(z0,x0,t)  -  A(z,x, t)e  c  (64) 

A(z,x,t)  is  the  envelope  of  the  pulse  and  is  the  carrier  freouency. 

The  nth  temporal  moment  is  defined  as 

CO 

/ 

Mn(z)  »  A*(z,x.t) tnA(z,x,t)dt  (65) 

—00 

where : 

n*0, 1 , 2 . . 
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4.2  Physical  Significance 


When  n*-0  Equation  (65)  becomes 


M0(z)  - 


|A(z,x,t) | 2dt 


(66) 


Equation  (66)  is  ju3t  the  expression  for  the  total  energy  in  the  pulse. 
When  n»l,  Equation  (65)  becomes 


Mi  (z) 


00 

(  A*(z,x,t)ta.(z,x,t)dt 


—00 


(67) 


Equation  (67)  is  the  expression  for  the  mean  of  the  square  of  the 
envelope.  If  Mj(z)  is  normalized  by  dividing  by  Mo(z),  then  an  expres¬ 
sion  for  the  mean  arrival  time  (denoted  by  t^)  is  obtained. 


ta  -  Mi(z)/M0(z) 


(68) 


The  mean  arrival  time  can  be  used  to  calculate  the  delay  caused  by  the 
electron  bubble  medium. 

When  moments  of  higher  order  are  needed,  t  is  used  to  define  the 
central  moment.  The  nth  central  moment  is  defined  as 


M  (z) 

n 


A*(z,x,t) (t-t  )  A(z,x,t)dt 


The  2nd  central  moment  is  analogous  to  the  variance  in  probability 
theory.  When  n*2.  Equation  (69)  (after  normalization)  becomes 


(69) 


M2(z)  M2(z) 

2  2 
M0(z)  “  1  "  M0(z)  _ta 


(70) 


t2  is  called  tne  mean  square  pulse  width,  x  gives  a  measure  of  the 
pulse  spreading  caused  by  frequency  dispersion  in  the  medium. 


The  third  central  moment  is  written  as 


M3(i)  -  j  A*(*<X*t)  (t-ta)  3A(Z,x,t)dt  (71) 

—00 

The*  skevness  s  is  defined  by 

M3(z) 

.  gT3  (72) 

M  O' (2) 

of 

st 3  *  H3{z)/M0^)-‘3M2(Ji)t  /(M0(i:))2+2t  3  (73) 

The  value  of  s  is  a  measure  of  the  pulse  asymmetry.  If  a  pulse  is 
perfectly  symmetric  about  the  mean  Arrival  time,  the  skewhess  Will  be 
zero.  In  this  thesis ^  if  s  is  negative,  more  of  the  energy  in  the 
puise  is  concentrated  in  the  trailing  edge  of  the  pulse  rather  in  its 
leading  edge» 


5. 


Da  ca 


In  this  section,  I  will  calculate  the  moments  of  the  pulse  envelope 
after  the  pulse  has  been  propagated  through  the  electron  bubble  medium. 
Pulses  of  three  different  carrier  frequencies  and  two  different  pulse 
widths  will  be  analyzed.  The  three  carrier  frequencies  are  2.5  GHz, 

800  MHz  and  500  MHz.  The  two  pulse  widths  are  6.51  nanoseconds  and 
12.95  nanoseconds.  For  each  of  these  carrier  frequencies  and  pulse 
widths,  a  receiver  will  be  placed  at  two  positions.  The  first  is 
x*5.297  km.  Note  from  Figure  2  that  at  x*5.297  km  the  bubble  has 
relatively  constant  electron  density.  The  second  x  position  is  x*8.856  km. 
At  this  x  position,  the  bubble  has  a  very  sharp  electron  density  fluc¬ 
tuation. 

SHEET  1  is  run  nine  times  for  each  carrier  frequency  at  each  x 
position  to  obtain  the  parabolic  equation  solution  across  the  entire 
pulse  bandwidth.  Linear  interpolations  are  done  to  obtain  the  ampli¬ 
tude  and  phase  of  the  parabolic  equation  solutions  between  the  SHEET  1 
outputs.  For  the  500  MHz  carrier  frequency,  additional  SHEET  1  runs 
were  made  because  the  amplitude  started  to  fluctuate  greatly  at  the 
lower  end  of  the  pulse  bandwidth.  Figures  5-16  show  the  interpolated 
parabolic  equation  solutions  for  each  carrier  frequency  and  each  x 
position.  These  Figures  are  essentially  the  transfer  function  which 
describes  the  electron  bubble  medium.  This  transfer  function  is  de¬ 
pendent  on  frequency  and  on  x. 

Some  of  the  input  and  output  pulses  are  shown  in  Figures  17-23. 

These  Figures  show  the  input  pulse,  the  input  pulse  spectrum,  the 
input  pulse  envelope,  the  output  pulse  at  x«5.297  km,  the  output  pulse 
envelope  at  x*5.297  km,  the  output  pulse  at  x»8.856  km  and  the  output 
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i^ure  5-  The  amplitude  of  the  transfer  function  for  2.5  GHZ 
at  x=5.297  km. 
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Figure  8-  The  phase  of  the  transfer  function  for  2.5  GHZ  at 
x=8.856  km. 


Figure  10-  The  phase  of  the  transfer  function  for  800  MHZ  at 
x-5.297  km. 
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Figure  11-  The  amplitude  of  the  transfer  function  for  800  MHZ 
at  x=8 . 856  km. 


Figure  12-  The  phase  of  the  transfer  function  for  800  MHZ  at 
x=8.856  km. 


Figure  15  The  amplitude  of  the  transfer  function  for  500  MHZ 
at  x=8.856  km. 


Figure  17-  The  input  pulse  applied  to  the  medium  with  a  c 
frequency  of  800  MHZ  and  a  pulse  width  of  6.51 
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Figure  19-  The  envelope  of  the  input  pulse. 


Figure  21-  The  output  pulse  envelope  at  x=5.297  km. 


pulse  envelops  for  x*8.856  km  respectively  for  the  6.51  nanosecond, 

800  MHz  pulse  (one  of  the  twelve  carrier  frequency  and  pulse  width 
combinations  analyzed).  In  Figure  17  (the  input  pulse),  the  leading 
edge  of  che  pulse  is  in  positive  time  and  the  trailing  edge  of  the 
pulse  is  in  negative  time.  The  pulse  is  centered  at  time  zero. 

Figures  22  and  23  clearly  show  an  example  of  pulse  distortion.  In 
Figure  22,  the  delayed  mean  arrival  time  (t  )  is  marked.  A  pulse  prop 
agating  through  free  space  the  same  distance  would  be  centered  at  t**Q. 
The  mean  arrival  time  is  the  excess  delay  caused  by  the  electron  bubbl 
medium. 

The  normalized  moments  of  the  output  pulse  envelopes  are  plotted 
versus  frequency  in  Figures  24-29.  Figures  24-26  are  for  the  pulse 
width  equal  to  6.51  nanoseconds  and  Figures  27-29  are  for  the  pulse 
width  equal  to  12.95  nanoseconds.  The  dashed  lines  are  the  data  for 
x-5.297  km  and  the  solid  lines  are  for  x*8.856  km. 
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In  section  2.1  the  fluctuating  part  of  the  relative  dielectric 


permittivity  was  derived  as  Equation  (29) 
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As  the  wave  frequency  is  made  larger,  e ; (r)  tends  toward  zero.  Thus 
for  high  frequencies  the  fluctuations  in  electron  density  will  have 
little  effect  on  a  propagating  wave.  The  wave  will  propagate  like  the 
medium  is  free  space.  However,  at  lower  frequencies,  the  electron 
density  fluctuations  will  have  a  great  effect.  Also  the  effect  will 
be  different  for  every  frequency  across  the  pulse  bandwidth.  Thus  a 
low  carrier  frequency  pulse  propagating  through  an  electron  bubble 
medium  will  be  distorted. 

To  analyze  the  results,  the  moments  of  Che  pulse  envelope  must 
be  related  to  its  frequency  spectrum.  The  moment  theorem  [Papoulis, 
<19b2)|,  relates  the  derivatives  of  F(w)  at  the  origin  to  the  moments 
of  the  time  function.  In  this  problem,  a  carrier  frequency  shifts 
the  frequency  spectrum  of  the  envelope  to  the  carrier  frequency  so 
the  moments  of  the  envelope  will  be  related  to  the  derivatives  of  its 
frequency  spectrum  at  the  carrier  frequency.  In  addition  to  the 
shifting,  the  frequency  spectrum  is  also  multiplied  by  a  factor  of  1/2. 

The  moment  theorem  which  relates  the  derivatives  of  the  frequency 
spectrum  to  the  moments  of  the  pulse  envelope  is 
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The  frequency  spectrum  of  the  envelope  is  written  in  terms  of  a  magni¬ 
tude  part  and  a  phase  part. 
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The  pulse  is  written  as  an  envelope  function  multiplied  by  a  complex. 
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Using  the  series  representation  for  the  exponential  gives 
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integrating  (83)  termwise  gives 
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and  .mbstituting  (74)  Into  (84)  results  in 
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Assuming  Mi)* l  and  equating  the  terms  of  (84)  and  (80)  results  in 
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and  defining 


tA(t)dt 


(89) 


gives 


<t-t  )2A(t)dt 
m 


(?0) 


(  :-t  )3A(t)dt 
tn 


(91) 


-db(u)  I 


u 


du 


(92) 


n2 


-d2H(w)  | 

diii^ 


(93) 


d  3d> (uj)  |  u“u) 

« --s? - £ 


(94) 


Thus  Che  delay  is  related  Co  Che  slope  u.f  the  phase  function  at  the 
carrier  frequency  and  the  pulse  spreading  is  related  to  the  curvature 
of  the  amplitude  function  at  the  carrier  frequency.  The  skewness  is 
related  to  the  3rd  derivative  of  the  phase  function.  Note  that  the 
delay,  pulsevldth  and  skc-w.iess  defined  above  are  slightly  different 
than  the  quantities  defined  in  Section  4.  The  definitions  above  are 
related  to  the  envelope  while  the  definitions  in  Section  4  are  related 
to  the  square  of  the  envelope  (since  the  pulses  are  real) .  In  the 
next  several  paragraphs,  the  above  derivations  will  be  used  to  qualita¬ 
tively  explain  the  output  pulses  due  to  several  different  transfer  func¬ 
tions.  The  delay,  pulsewidth  and  skewness  of  the  square  of  the  envelope 
are  directly  related  to  the  delay,  pulsewidth  and  skewness  of  the  envelope 
itself . 
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According  to  Equation  (29),  one  would  expect  pulse  distortion  to 
become  greater  when  the  carrier  frequency  is  decreased.  Thus  the  pulse 
width  and  the  skewness  should  becoae  greater  when  the  carrier  frequency 
Is  decreased.  That  is  exactly  what  is  observed  in  Figures  25,  26,  28, 
and  29  when  the  receiver  is  placed  at  x-5.297  km.  But  when  the  re¬ 
ceiver  is  placed  at  x»8.856  km,  something  different  is  observed.  From 
2.5  GH2  to  800  MHz  the  expected  distortion  versus  frequency  occurs.  But 
at  500  MHz  the  distortion  becomes  less.  By  studying  Figures  12  and  16 
one  can  see  why  the  distortion  is  less  at  the  lower  frequency.  In 
Figure  12,  the  800  Mhz  phase  transfer  function,  there  is  a  non-linearity 
near  the  center  of  the  bandwidth.  Xn  Figure  16,  the  500  MHz  phase 
transfer  function,  the  phase  is  nearly  linear  near  the  center  of  the 
pulse  bandwidth.  According  to  Equation  (92),  non-linear  phase  causes 
pulse  distortion  because  delay  is  dependent  on  the  derivative  of  the 
phase  function.  Delay  as  a  function  of  frequency  (D(f)  is  written  as 

0(f)  -  (95) 

If  the  phase  is  non-linear,  some  of  the  frequency  components  in  the 
pulse  are  delayed  different  amounts  of  time  than  other  components. 

Thus  the  pulse  is  distorted.  So  one  reason  the  800  MHz  pulse  is  dis¬ 
torted  more  than  the  500  MHz  pulse  is  because  the  phase  is  more  non¬ 
linear  near  the  center  of  the  800  MHz  phase  transfer  function. 

Figures  11  and  15  can  be  used  to  explain  why  the  pulsewidth 
decreases  at  500  MHz  for  x"8,856  km.  When  the  mrmUrjde  spectrum 
of  the  input  pulse  is  multiplied  by  the  transfer  function  of  Figure  1). 
a  large  amount  of  negative  curvature  at  the  carrier  frequency  results. 
However  when  the  amplitude  spectrum  of  the  input  pulse  is  multiplied 
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by  the  transfer  function  of  Figure  15,  less  negative  curvature  at 
the  carrier  frequency  results  because  in  Figure  15  there  is  a  large 
amount  of  positive  curvature  which  tends  to  cancel  some  of  the  negative 
curvature  of  the  input  spectrum.  According  to  Equation  (93) ,  the 
pulsewidth  is  related  to  the  negative  curvature  of  the  amplitude 
function  at  the  carrier  frequency.  Thus  the  output  pulse  due  to  the 
transfer  function  of  Figure  11  will  be  spread  in  time  more  than  the 
output  pulse  due  to  the  transfer  function  of  Figure  15. 

Figure  12  can  also  be  used  to  justify  the  large  skewness  value 
at  x"8.856  km  for  the  800  MHz  carrier  frequency.  Since  skewness  is 
related  to  the  3rd  derivative  of  the  phase  function  according  to 
Equation  (94),  the  transfer  function  of  Figure  12  would  give  rise 
to  a  large  skewness  value  because  of  the  non-linearities  in  the  phase 
near  the  center  of  the  bandwidth. 

At  first  the  above  results  might  be  considered  disturbing  because 
Equation  (29)  predicts  that  the  pulse  distortion  should  ri3e  as  the 
carrier  frequency  is  decreased.  However,  more  analysis  could  easily 
rationalize  the  results.  If  a  receiver  was  placed  at  every  x  point 
below  the  bubble  and  it  a  pulse  was  applied  at  every  x  point  above 
the  bubble,  on  tne  average  the  expected  pulse  distortion  would  occur. 

In  the  real  ionosphere,  the  electron  bubble  would  be  moving 
ijlative  to  the  line  of  propagation.  If  a  pulse  train  was  transmitted, 
the  delay,  pulsewidth  and  the  skewness  of  each  successive  pulse  would 
vary.  But  on  the  average,  the  delay,  pulsewidth  and  skevnesa  would 
rise  if  the  carrier  frequency  of  the  pulse  train  was  decreased. 

In  order  to  use  this  computer  model  to  study  this  propagation 
problem  in  greater  detail,  a  receiver  would  have  to  be  placed  at 
every  x  position.  To  do  that  would  require  a  very  large  computer  budget. 
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The  software  procedure  used  to  obtain  the  moments  for  a  pulse 
envelope  is  shown  below.  Each  of  tha  programs  is  explained  and  listed 
in  this  appendix. 

1.  Program  SHEET  1 

This  program  numerically  solves  the  parabolic  equation  for  a 
wave  propagating  through  the  model  electron  bubble.  Program  SHEET  1 
is  explained  fully  in  section  3.2. 


2.  Program  Pulse 

This  program  computes  the  frequency  spectrum  of  the  Incident 
pulse  from  the  time  function.  The  frequency  spectrum  is  computed 
using  the  Discrete  Fourier  Transform  via  a  standard  decimation  in 
time  FFT  routine. 

The  pulse  time  function  is  set  up  In  an  array.  The  frequency 

spectrum  is  computed  using  the  DFT. 

_  -*(  >mn 
F. (mo)  •  l  p  (nT)e 
n-0  1 

where 


N  -  8192 

Pi(nT)  -  the  pulse  time  function 
m  •  0,1,2, . . .N-l 
T  -  sampling  rate 

The  frequency  spectrum  is  computed  for  N  values  in  the  range  mp»(0,2ir) 
radians.  To  obtain  the  actual  frequency  scale,  the  following  relation- 


So 


F^(tno)  •  F^(majT) 

In  this  thesis,  I  will  adjust  T  each  time  the  carrier  frequency  is 
changed  so  thac  the  1024th  point  of  the  frequency  spectrum  array  cor¬ 
responds  to  the  carrier  frequency. 

3.  Program  INTER 

This  program  interpolates  between  the  nine  SHEET  1  solutions  to 
obtain  a  parabolic  equation  solution  corresponding  to  each  of  the 
points  in  the  frequency  spectrum  computed  by  program  PULSE.  Two  linear 
interpolations  are  done,  one  for  the  amplitude  and  one  for  the  phase. 

4.  Program  MULT 

This  program  multiplies  the  interpolated  parabolic  equation  solu¬ 
tions  and  the  input  pulse  frequency  spectrum  to  obtain  the  output 
pulse  frequency  spectrum. 

5.  Program  IFFT 

Program  IFFT  calculates  the  output  pulse  time  function  from  the 
frequency  spectrum  using  the  inverse  Discrete  Fourier  Transform  via 
a  standard  decimation  in  time  FFT  routine.  The  inverse  DFT  is  given 
by 

.  N-l  id~)n.n 

Po(nT)  -  jj-  >  F0(mo)e 
m-0 

where : 

N  -  8192 

n  ■  0,1,2,.... N-l 

Fg(mo)  *  output  pulse  frequency  spectrum 


Id  7 


Po(nT)  ■  output  pulse  time  function 
T  ■  sampling  rate 


6.  Program  ENVEL 

Pregram  ENVEL  finds  the  envelope  of  the  output  pulse  and  then 
calculates  the  first  four  moments. 
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PROGRAM  SHEETl ( INPUT , OUTPUT t ELFL , RESULT , PASS , DA2 , 

+TAPE1 “RESULT , TAPE2  -ELFL , TAPE3  -DA2 , TAPE5-PASS , 

+TAPE6 -OUTPUT) 

****************************************************************** 

..THIS  PROGRAM  SOLVES  THE  PARABOLIC  EQUATION  NUMERICALLY 
..FOR  A  WAVE  PROPAGATING  IN  THE  IRREGULAR  IONOSPHERE 
..OR  IN  FREE  SPACE. 

**************•**#★***★*****★*★★*★******★★**#*<,'*****#************** 


..THE  4  INPUT/OUTPUT  FILES  PREFORM  THE  FOLLOWING  FUNCTIONS: 

(1)  DA2  CONTAINS  INPUT  PARAMETERS 
K,N,TAU,H,SIGM,FREQ,KSKIP, 

XM,KK,KL,AL,BET, 20 ,  ZOO  ,  EMAX ,  EZ  ,  ZMAX . 

(2)  ELFL  CONTAINS  INFORMATION  ABOUT  THE  HORIZONTAL 
ELECTRON  DENSITY  STRUCTURE. 

(3)  RESULT  CONTAINS  THE  SOLUTIONS  TO  THE  PARABOLIC 
EQUATION  AT  EACH  HEIGHT  Z  AND  EACH  X  STEP. 

(4)  PASS  IS  THE  RENAMED  FILE  RESULT.  PASS  IS  USED 
IF  FURTHER  COMPUTATIONS  ARE  TO  BE  DONE  IN  FREE 
SPACE  BELOW  THE  BUBBLE. 

****************************************************************** 

★★★★★a************************************************************ 


...THE  4  SUBROUTINES  PREFORM  THE  FOLLOWING  FUNCTIONS 

(1)  UO  SETS  THE  INTIAL  VALUES  ON  THE  WAVE  AT  Z-0. 

(2)  FLUCT  CALCULATES  THE  FLUCTUATING  PART  OF  THE 
DIELECTRIC  PERMITTIVITY  AND  CALCULATES  THE  F 
FUNCTION. 

(3)  DCALC  CALCULATES  THE  MATRIX  EQUATION  D-A*U. 

(4)  THOMAS  SOLVES  THE  MATRIX  EQUATION  B*X-D  USING 
GAUSSIAN  ELIMINATION  VIA  THE  THOMAS  ALGORITHM. 

♦★★♦★a************************************************************* 

...THE  INPUT  PARAMETERS  IN  DA2  ARE: 

(1)  K  THE  NUMBER  OF  STEPS  IN  THE  X  DIRECTION. 

(2)  N  THE  NUMBER  OF  STEPS  IN  THE  Z  DIRECTION. 

(3)  H  THE  STEP  SIZE  IN  THE  X  DIRECTION  IN  METERS. 

(4)  TAU  THE  STEP  SIZE  IN  THE  Z  DIRECTION  IN 

METERS . 

(5)  SIGM  A  FACTOR  USED  IN  THE  NUMERICAL  PARABOLIC 
EQUATION  SOLUTION.  (SEE  REPORT) 

(6)  FREQ  THE  WAVE  FREQUENCY  IN  MHZ. 

(7)  KSKIP  THE  NUMBER  OF  STEPS  IN  THE  Z  DIRECTION 
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RECORDED  ON  FILE  PASS.  IF  KSKIP  IS  NOT  EQUAL 
TO  ZERO  THEN  PASS  IS  AUTOMATICALLY  READ  AND 
ELFL  IS  NOT  READ. 

(8)  3M  ONLY  EVERY  ZM  STEP  IN  THE  Z  DIRECTION  IS 
RECORDED  ON  FILE  RESULT. 

(9)  KK  ONLY  EVERY  KK  STEP  IN  THE  X  DIRECTION  IS 
RECORDED  ON  FILE  RESULT. 

(10)  KL  THE  NUMBER  OF  STEPS  IN  THE  Z  DIRECTION 
COMPUTED  IN  ALL  PREVIOUS  RUNS. 

(11)  AL  A  FACTOR  USED  IN  THE  WEIGHTING  FUNCTION. 

(12)  BET  A  FACTOR  USED  IN  THE  WEIGHTING  FUNCTION. 

(13)  ZO  A  FACTOR  USED  IN  THE  WEIGHTING  FUNCTION. 

(14)  ZOO  A  FACTOR  USED  IN  THE  WEIGHTING  FUNCTION. 

(15)  EMAX  A  FACTOR  USED  TO  COMPUTE  THE  PARABOLIC 
PROFILE.  THE  ELECTRON  DENSITY  IN 
ELECTRONS/CM* *3  AT  ZMAX. 

(16)  EZ  THE  ELECTRON  DENSITY  IN  ELECTRONS/CM* *3 
AT  7,-0. 

(17)  ZMAX  THE  HEIGHT  OF  THE  MAXIMUM  ELECTRON 
DENSITY  IN  THE  PARABOLIC  PROFILE. 

******************************************************************* 

******************************************************************* 

COMMON /AAA/ELX ( 800 ) 

COMMON/ BBB/FREQ ,H,N, TAU , PI , XKI 
COMMON/CCC/AL , BET , EMAX , EZ , ZMAX , ZO , ZO  0 
COMPLEX  AK ( 8 0 0 )  ,CK(80O)  # U ( 800 )  ,CCK(800)  ,D(800) 

+ , X ( 800 ) , AAK (800) 

DIMENSION  F ( 800 )  ,AMP(800)  ,PH(800)  ,TEMP(2,800) 

INTEGER  ZM 

. . . READ  THE  INPUT  PARAMETER  FROM  FILE  DA2 . 

READ (3, 100) K,N, TAU, H,SIGM, FREQ, KSKIP 
READ (3,125) ZM,KK,KL 

READ ( 3 , 150 ) AL , BET, Z0 , Z00 , EMAX , EZ , ZMAX 

...NORMALIZE  THE  STEP  SIZE  TO  THE  WAVE  NUMBER  AND  SET  THE 
...VALUE  FOR  PI.  THIS  IS  EQUATION  (35). 

PI-4.0  *ATAN ( 1 . ) 

XKI-300. / ( 2*PI*FREQ) 

TAU-TAU/XKI 

H-H/XKI 

...IF  KSKIP  IS  NOT  EQUAL  TO  ZERO  THEN  THE  COMPUTATIONS  WILL 
...BE  DONE  OUTSIDE  THE  IRREQULAR  MEDIUM.  THE  FILE  ELFL 
...DOES  NOT  NEED  TO  BE  READ  AND  THE  INITIAL  VALUES  FOR 
...THE  l"  MATRIX  ARE  OBTAINED  FROM  THE  PREVIOUS 
...COMPUTATIONS  IN  FILE  PASS. 

IF { KSKIP . NE. 0 ) GO  TO  15 
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C...READ  THE  FILE  ELFL  WHICH  CONTAINS  THE  ELECTRON  DENSITY 
C ...  INFORMATION  OF  FIGURE  2.  THIS  INFORMATION  IS  STORED 
C. . .IN  ELX. 

C 

READ (2  , 175)  (ELX(I)  , I«1,K) 

...IF  FILE  PASS  HAL  NO  VALUES  STORED  IN  IT,  THE  COMPUTATIONS 
...START  AT  Z«0  AND  THE  PARAMETERS  ARE  OUTPUT  TO 
...FILE  RESULT. 

KL-N 

T*TAU*XKI 
Hl-H*XKI 

WRITE (1, 200) K,N,T,H1,SIGM, FREQ, KL,KSKIP 
WRITE (1,225)  KK,ZM 

...THE  U  MATRIX  IS  THE  COMPLEX  AMPLITUDE  OF  THE  WAVE. 
...SUBROUTINE  UO  IS  CALLED  TO  SET  THE  INITIAL  CONDITIONS 
. . .OF  U. 

CALL  UO ( U , K ) 

GO  TO  11 
C 

C. . .THIS  IS  THE  BRANCH  POINT  IF  KSKIP  IS  NOT  EQUAL  TO  ZERO 
C...AND  THE  FILE  PASS  NEEDS  TO  BE  READ.  IN  THIS  CASE  THE 
C...INTIAL  VALUES  OF  U  ARE  OBTAINED  FROM  PREVIOUS 
C.  .  ..COMPUTATIONS  STORED  IN  FILE  PASS. 

C 

15  READ (5,250)  K , N, TAU, H, SIGM, FREQ , KL, KSKIP 
READ (5,275)  KK, ZM 

DO  16  L«l, KSKIP 

READ (5, 300) (TEMP (1,1) , TEMP (2, I) ,AMP(I) ,PH(I) ,I«1,K) 

16  CONTINUE 

DO  17  I«1,K 

U ( I ) »CMPLX ( TEMP (1,1) , TEMP (2,1) ) 

17  CONTINUE 
11  CONTINUE 

C 

C...THE  NUMERICAL  SOLUTION  TO  THE  PARABOLIC  EQUATION  BEGINS. 

C 

C...THE  TWO  MATRICES  OF  EQUATION  (43)  A  AND  B  ARE  STORED  IN 
C... THREE  ARRAYS  EACH.  FOR  A,  THE  THREE  ARRAYS  ARE  AK , CK , BK . 
C...AK  AND  BK  STORE  THE  OFF  DIAGONAL  ELEMENTS  OF  A 
C...AND  CK  STORES  THE  DIAGONAL  ELEMENTS  OF  A.  SIMILARILY, 
C...AAK,  BBK ,  CCK  STORE  THE  ELEMENTS  OF  MATRIX  B. 

C... SINCE  THIS  PROBLEM  HAS  SPECIAL  SYMMETRY  AK=BK  AND 
C . . . AAK*BBK  SO  THESE  ARRAYS  ARE  NOT  STORED  TWICE. 
r 

C...THE  AK  AND  AAK  MATRICES  ARE  SET  UP  ACCORDING  TO  EQUATIONS 

C.  .  .  (45-48)  . 
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DO  1  1*1  /  K 
R-0.0 

AI*0 ♦ 5/ (H*H) 

AK  1 1)  -CMPLX(R,AI) 

R*0 . 0 

AI* ( . 5/ (H*H) ) *TAU*SIGM 
AAK(I)  “CMP  LX  ( R,  AI) 

1  CONTINUE 

C 

C...NOW  THE  MAIN  LOOP  BEGINS.  FIRST  THE  DIAGONALS  OF  THE 
C...A  AND  B  MATRICES  ARE  SET  UP.  BY  EQUATIONS  (45-48) 

C ...  IT  CAN  BE  SEEN  THAT  THESE  VALUES  DEPEND  ON  THE  ELECTRON 
C... DENSITY  FLUCTUATIONS.  IF  THE  COMPUTATIONS  ARE  TO  BE 
C . . . DONE  IN  FREE  SPACE  THERE  ARE  NO  ELECTRON  DENSITY 
C ...  FLUCTUATIONS  SO  THE  F.  FUNCTION  IS  SET  TO  ZERO. 

C 

DO  2  J*1,N 

DO  3  1*1 ,K 
F ( I) *0.0 
3  CONTINUE 

IF ( KSKIP . NE. 0 )  GOTO  5 
CALL  FLUCT(F,J,K) 

5  DO  6  1*1 , K 
R-0.0 

AI— 1.0/(H*H)  +F(I)/2 
CK  ( I)  ■  CMP  LX  ( R,  AI) 

R*l .  0 

AI—  (TAU*SIGM)  /(H*H)  +F  ( I)  *SIGM*TAU/2 
CCK(I) “CMPLX (R, AI) 

6  CONTINUE 

...NEXT  THE  PRODUCT  OF  TWO  MATRICES  D-A*U  IS  CALCULATED  IN 
... DCALC . 


CALL  DCALC ( AK  /  AK ,  CK ,  U ,  D ,  K ) 

...NOW  THE  SYSTEM  OF  EQUATIONS  IS  SOLVED  VIA  THE  THOMAS 
...ALGORITHM.  X-(B  INVERSE) *A*U.  THIS  IS  EQUATION  (51). 

CALL  THOMAS(AAK, AAK , CCK, X , D , K) 

...FINALLY  THE  EQUATION  U (N+l ) -U (N) -X*TAU  IS  CALCULATED. 
...EQUATION  (49). 

DO  4  1*1 , K 

R*R£AL(U(I) ) -TAU*REAL (X (I) ) 

AI-AIMAG(U(I) ) -TAU* AIMAG ( X ( I ) ) 

U(I)  -CMPLX  (R,AI) 

4  CONTINUE 

C 

C...THE  MAGNITUDE  AND  PHASE  OF  THE  COMPLEX  AMPLITUDE  IS 
C.  . .CALCULATED. 


o  o  o  n  o  o 


DO  7  I-1,K 

AMP ( I) -SQRT ( REAL ( U  ( I ) ) **2+AIMAG (U { I)  ) **2) 
PH( I) -ATAN2 (AIMAG (U ( I) ) , REAL ( U ( I ) ) ) 

I P { PH { I ) . LT . 0 . 0 )  PH (I) -PH (I) +2*PI 
PH  ( I)  -PH  ( I)  /PI 
7  CONTINUE 

C 

C... WRITE  THE  REAL  PART/  THE  IMAGINARY  PART,  THE  AMPLITUDE 
C...AND  THE  PHASE  TO  FILE  RESULT.  ONLY  EVERY  ZM  TH  STEP 
C ...  IN  THE  Z  DIRECTION  AND  EVERY  RK  TH  STEP  IN  THE  X 
C. . .DIRECTION  IS  OUTPUT, 

C 

I- ( J/ZM) *ZM 
IF(I.NE.J)  GOTO  2 

WRITE ( 1 , 300) (REAL(U(I) ) ,AIMAG(U(I) ) ,AMP(I) 

♦  , PH ( I) , 1-1 , K , KK) 

2  CONTINUE 

ENDFILE  1 

100  FORMAT ( 14 , lX, 14,4 (1X,F12.6) ,1X,I6) 

125  FORMAT (14,14,16) 

150  FORMAT (7E10. 2) 

175  FORMAT (10F8. 2) 

200  FORMAT ( 2 HR- , 14 , 2X , 2HN- , 15 , 2X , 4HTAU- , FI 1 . 4 , 2X 
+  , 2HH-, Fll. 4 , 2X, 5HSIGM-, F6 . 4 , 2X , 5HFREQ-, F8 . 2 , 2X 
+  ,  3HRL-,  15 , 2X,  6HKSRIP-,  IS) 

225  FORMAT (3HRK-, 14, 3HZM-, 14) 

250  FORMAT (2X,I4,4X,I5/6X,F11.4,4X/F11.4,7X,F6.4,7X,F8.2, 
+5X,IS,8X,I5) 

275  FORMAT (3X, 14, 3X, 14) 

300  FORMAT (2 (2X, 4E16 . 10) ) 

400  FORMAT ( IX, 15 ) 

STOP 
END 


SUBROUTINE  UO(U,L) 

...THIS  SUBROUTINE  SETS  THE  INITIAL  VALUES  TO  THE  MATRIX  U 

COMPLEX  U ( L) 

DC  1  I-1,L 

UC)  -CMPLX  (1.0/C.0) 

l  cownruE 

RETURN 

END 

C 

£★*★************************************************************** 

C 

SUBROUTINE  FLUCT(F,J,K) 

C 


1  ‘>3 


C . . . THIS  SUBROUTINE  CALCULATES  THE  FLUCTUATING  PART  OF  THE 
C. . .DIELECTRIC  PERMITTIVITY  AND  THE  F  FUNCTION  FOR  THE  BUBBLE 
C... MODEL.  THE  VERTICAL  VARIATIONS  ARE  DETERMINED  BY  THE 
C. . .PARABOLIC  PROFILE  MULTIPLIED  BY  THE  WEIGHTING  FUNCTION 
C...G.  THE  HORIZONTAL  VARIATIONS  ARE  DETERMINED  BY  THE  FILE 
C...ELFL  (STORED  IN  ARRAY  ELX) . 

COMMON /AAA/ELX (800) 

COMMON /BBB/FREQ , H , N , TAU , PI , XKI 
COMMON/CCC/AL , BET , EMAX , EZ , ZMAX , Z0 , Z0  0 
DIMENSION  F(600) ,DPZ(600) ,DPZ1 (600) 

Z«TAU*(J-1) *XKI 
Zl«TAU*J*XKI 
C 

C. . .CALCULATE  THE  WEIGHTING  FUNCTION  FOR  Z  AND  Zl.  (G  AND  Gl) 
C 

IF(Z.LE.ZO)  GOTO  1 
IF(Z.GT.ZOO)  GOTO  2 
G*1 . 0 
Gl-1.0 
GO  TO  3 

1  ARG»( (Z-ZO) /AL) **2 
ARGl-( (Z1-Z0) /AL) **2 
IF(ARG.GT. 200.0)  GOTO  4 
IF (ARG1.GT. 200.0)  GOTO  4 
G- EXP(-ARG) 

Gl-EXP(-ARGl) 

GOTO  3 

2  ARG* ( ( Z-Z00) /BET) **2 
ARGl-( (Zl-ZOO) /BET) **2 
IF(ARG.GT. 200.0)  GOTO  4 
IF(ARGI.GT. 200.0)  GOTO  4 
G-EXP(-ARG) 

Gl-EXP(-ARGl) 

GOTO  3 
4  G-0 . 0 

Gl-0.0 
C 

C...THE  SCALE  HEIGHT  HH  IS  CALCULATED. 

C 

3  HH-110.0E+03 
C 

C...THE  ELECTRON  DENSITY  AT  THE  EDGE  OF  THE  BUBBLE  IS 
C ...  CALCULATED  FOR  HEIGHTS  Z  AND  Zl. 

C 

E«EMAX* (1.0-( (Z-ZMAX) /HH) **2) *G 
El-EMAX* (1. 0- ( (Zl-ZMAX) /HH) **2) *G1 
FACT- I . 0 
FACT1-1 . 0 

IF(E.LE.O.O)  FACT-0.0 
IF (El. LE. 0.0)  FACT1-0 . 0 
C 

C...THE  SQUARE  OF  THE  PLASMA  FREQUENCY  AND  THE  FACTOR 


I1. 4 


C. . . (WP/W) **2  ARE  CALCULATED  FOR  THE  BUBBLE  EDGE  AT 
C... HEIGHTS  Z  AND  Zl. 

C 

FP-80.6E-6*E/(FREQ**2) 

FP1»80.6E-6*E1/(FREQ**2) 

C 

C...THE  FLUCTUATING  PART  OF  THE  DIELECTRIC  PERMITTIVITY 
C...AND  THE  F  FUNCTION  IS  CALCULATED  FOR  EACH  X  STEP  AT 
C.. .HEIGHT  Zl.  EQUATION  (32). 

C 

DO  6  1*1 , K 

DPZ(I) -FACT*FP* (ELX ( I) /ELX (K) -1.0) 

DPZ1 (I) -FACT1*FP1* (ELX (I) /ELX (K) -1.0) 

F  ( I)  -  (DPZ  ( I)  +DPZ1  ( I)  )/2.0 
6  CONTINUE 

RETURN 
END 
C 

C* **************************************************************** 

C 

SUBROUTINE  DCALC ( AK , BK , CK , U , D , K ) 

C 

C. . .THIS  SUBROUTINE  CALCULATES  THE  PRODUCT  D-A*U.  EQUATION 
C...(54).  NOTE  THAT  AK-BK  IN  THIS  PROBLEM. 

C 

COMPLEX  AK(K)  ,BK(K)  ,CK(K)  ,U(K)  ,D(K) 

M-K-l 

DO  1  1-2, M 

D  ( I)  -AK  ( I)  *U(I-1)  +CK(I)  *U(I)  +BK(I)  *U(I+1) 

1  CONTINUE 

D(1)-CK(1)  *U(1)  +BK  ( 1)  *U  (2) 

D  (K)  -AK  (K)  *U  (K-l )  -t-CK(K)  *U(K) 

RETURN 

END 

C 

c* ************************** ************************************** 

C 

SUBROUTINE  THOMAS (A, B , C , X , D, K) 

COMPLEX  A ( 600 ) ,B(600) ,C<600) ,X(600) ,D(600) 

+ , ALF ( 600 ) , BT ( 600 ) 

C 

C . . . THIS  SUBROUTINE  SOLVES  THE  MATRIX  EQUATION  BX-AU 
C... (WHICH  IS  THE  SAME  AS  BX-D)  FOR  X. 

C 

C... CONSIDER  THE  SYSTEM: 

C 

C  (1)  A(I)  *X(I-1) +C(I)  *X(I)+B(I)*X(I  +  1)-D(I) 

C  (2)  X ( 1 ) — B ( 1 ) /C ( 1 ) *X ( 2 ) +D ( 1 ) /C ( 1 ) 

C  (3)  X (K) — A(K) /C(K) *X(K-1) +D(K) /C(K) 

C 

C... BECAUSE  X ( 0 ) *0  AND  X(K+l)-0 


ooooooooooo 


...ASSUME  THE  SOLUTION: 

X(I)-ALP(I)  *X(I+1)  +BT(I+1) 


. . . WHERE : 

(4)  ALF (2) »-B ( 1) /C (1)  SEE  EQ.  (2)  ABOVE 

(5)  BT(2)-D(1)/C(1)  SEE  EQ.  (2)  ABOVE 

(6)  ALFU+l)  »-B  (I)  / (A ( I)  *ALF(I)  +C(I)  ) 

(7)  BT  ( 1+1)  *  (D  ( I)  -A  (I)  *BT  (I)  )  /  (A{  I)  *ALF(I)  +C(I)  ) 

...NOTE  THAT  A»B  IN  THIS  PROBLEM  BECAUSE  OF  SYMMETRY. 

ALF (2) *~B (1) /C ( 1) 

BT(2)*D(1)/C(1> 

M*K-1 

DO  1  1-2, M 

ALF  (Id)  *-B  (I)  /  ( A  ( I)  *ALF(I)  +C(I)  ) 
BT(I+1)-(D(I)-A(I)  *BT (I)  )  / ( A(  I)  *ALF (I)  +C  ( I)  ) 

1  CONTINUE 

X  (K)  *  (D  (K)  -A (K)  *BT (K)  )  /  (A(K)  *ALF(K)  4-C  ( K }  ) 

DO  2  I»1,M 
J-M-I+l 

X ( J) “ALF { J+l) *X ( J+I) +BT( J+l) 

2  CONTINUE 
RETURN 

END 


i  > 


oono  oon  ononrj  ooo  oooo  noonno 


PROGRAM  PULSE ( INPUT , OUTPUT , I SPECT , TAPE 2 0 - 1 S PECT , I FREQ , 
+TAPE30*IFREQ) 

...THIS  PROGRAM  WILL  CALCULATE  AND  PLOT  THE  FFT  OF  A  PULSE 
.  . . THE  PARAMETERS  OF  THE  PULSE  THAT  CAN  BE 
.  ..CHOOSEN  ARE: 

1)  FREQ,  THE  CARRIER  FREQUENCY  IN  GHZ. 

2)  TAU,  THE  PULSE  DURATION  IN  NANOSEC. 

3)  A,  THE  AMPLITUDE  OF  THE  PULSES  IN 
VOLTS. 

...THE  PULSE  TIME  FUNCTION  IS  ALSO  PLOTTED. 

...THE  MAX  FREQUENCY  COMPONENT  SHOULD  BE  LESS  THAN  60  GHZ 
...TO  AVOID  ALIASING. 

DIMENSION  DATA (8200) ,XT(8200) 

COMPLEX  C(8200) 

INTEGER  PERIOD,  TAU 

...SET  UP  THE  PARAMETERS  OF  THE  PULSE. 

READ (30, 125)  FREQ, TAU 
A*l.  0 

...INITIALIZE  THE  ARRAYS,  DATA  (WHERE  THE  TIME  FUNCTION  AND 
...LATER  THE  FREQ  SPECTRUM  IS  STORED).  C  THE  COMPLEX  REP. 
...OF  DATA.  AND,  XT  AN  ARRAY  USED  TO  PLOT  AGAINST. 

DO  1  1*1,8200 
C (I) -CMPLX (0.00, 0.00) 

DATA (I) *0.0 
XA-I 

XT  (I)  *  (XA/  (FREQ*  8 .0)-4096.0/  (FREQ  *8 . 0 )  ) 

1  CONTINUE 

...NOW  THE  PULSE  IS  COMPUTED  AND  STORED  IN  DATA  AND  C. 

TAU1"TAU*(8*FREQ) 

J* (4096-TAU1/2) 

L* (4096+TAU1/2) 

DO  3  J1*J, L 
X-(Jl-J) *(3.14)/4.0 
DATA(Jl) *A*COS (X) 

3  CONTINUE 

...THE  PULSE  IS  MULTIPLIED  BY  A  GAUSSIAN  FUNCTION'  SO  THAT  IT 
...DOESN'T  BEGIN  AND  END  ABRUPTLY. 

IM*0 

V*SQRT( ( (TAU/2 ) **2) /(2*4 , 61) ) 

DO  8  I-J,L 

X*( (X T(I) -IM) **2) /(2*(V**2)  ) 

DATA ( I) *DATA ( I ) *EXP(-X) 

8  CONTINUE 


IH7 


DO  10  1-1,8192 
Yl-DATA(I) 

Y2-0 .0 

C(I) -CMPLX ( Yl , Y2 ) 

10  CONTINUE 

C 

C...PIND  THE  POWER  IN  THE  TIME  DOMAIN  EXPRESSION  OF  THE 
C. . .PULSE. 

C 

PT-0 . 0 

DO  15  I-J , L 
PT-DATA ( I ) * *2+PT 
15  CONTINUE 

C 

C . . . PLOT  THE  TIME  FUNCTION. 

C 

Xl-TAU/2+30 
CALL  USTART 
CALL  UERASE 
CALL  URESET 

CALL  UDAREA (0.00,7.49,0.00,5.71) 

CALL  USET ( "XBOTH" ) 

CALL  USET ( " YBOTH" ) 

CALL  UPSET ("XLABEL", "TIME  IN  NANOSEC.;") 

CALL  UPSET (“YLABEL", "AMPLITUDE;") 

CALL  USET ( "OWNSCALE" ) 

CALL  UWINDO (-Xl, XI, -1.0, 1.0) 

CALL  UPLOT1 (XT, DATA, 8192.0) 

CALL  UPAUSE 
C 

C...NOW  THE  FFT  SUBROUTINE  IS  CALLED,  THE  FREQ  SPECTRUM  WILL 
C ...  BE  RETURNED  IN  ARRAY  C. 

C 

9  NN-13 

CALL  FFT(C,NN, .TRUE. ) 

C 

C. . .NEXT  THE  FREQUENCY  SPECTRUM  IS  PLOTTED. 

C 

Xl-FREQ/1024 

DO  4  1-1,4097 

DATA ( I ) -SQRT { ( REAL ( C ( I ) ) **2) + (AIMAG (C ( I) ) **2) ) 

XT ( I) • ( 1-1) *X1 
4  CONTINUE 

C 

C. . . FIND  THE  POWER  IN  THE  FREQUENCY  DOMAIN  EXPRESSION  OF  THE 
C... PULSE  AND  ADJUST  THE  FREQUENCY  SPECTRUM  SO  THAT 
C. . .PARSEVAL'S  THEOREM  HOLDS. 

C 

PF-0.0 

DO  14  1-1,4096 

DATA (I) -DATA ( I) /I 500.0 

PF-DATA ( I ) **2+PF 


DO  10  1-1,8192 
Yl-DATA(I) 

Y2-0 .0 

C(I)  "CMPLX  ( Yl , Y2) 

10  CONTINUE 

C 

C. . . FIND  THE  POWER  IN  THE  TIME  DOMAIN  EXPRESSION  OF  THE 
C. . .PULSE. 

C 

PT-0 . 0 

DO  15  I«J,L 
PT-DATA ( I ) **2+PT 
15  CONTINUE 

C 

C . . . PLOT  THE  TIME  FUNCTION. 

C 

Xl-TAU/2+30 
CALL  US TART 
CALL  UERASE 
CALL  URESET 

CALL  UDARE A (0.00,7.49,0.00,5.71) 

CALL  USET ( "XBOTH" ) 

CALL  USET("YBOTH") 

CALL  UPSET ("XLABEL", "TIME  IN  NANOSEC.;") 

CALL  UPSET ("YL ABEL", "AMPLITUDE;") 

CALL  USET ( " OWN SCALE " ) 

CALL  UWINDO ( -XI , Xl , -1 . 0 , 1 . 0 ) 

CALL  UPLOT1 (XT, DATA, 8192,0) 

CALL  UPAUSE 
C 

C...NOW  THE  FFT  SUBROUTINE  IS  CALLED,  THE  FREQ  SPECTRUM  WILL 
C ...  BE  RETURNED  IN  ARRAY  C. 

C 

9  NN-13 

CALL  FFT (C,NN, .TRUE. ) 

C 

C . . . NEXT  THE  FREQUENCY  SPECTRUM  IS  PLOTTED. 

C 

Xl-FREQ/1024 

DO  4  1-1,4097 

DATA ( I ) -SQRT { (REAL(C(I) ) **2) + ( AIMAG (C ( I) ) **2) ) 
XT(I)-(I-1)  *X1 
4  CONTINUE 

C 

C. . . FIND  THE  POWER  IN  THE  FREQUENCY  DOMAIN  EXPRESSION  OF  THE 
C... PULSE  AND  ADJUST  THE  FREQUENCY  SPECTRUM  SO  THAT 
C. . .PARSEVAL'S  THEOREM  HOLDS. 

V- 

PF-0.0 

DO  14  1-1,4096 

DATA ( I ) -DATA (I)  /1 500.0 

PF-DATA ( I ) **2+PF 


1  '■'< 


14  CONTINUE 

XF-SQRT (PT/PF) 
XF-XF/2.5066 

DO  13  1*1 , 4096 
DATA (I) -DATA (I) *X? 


13  CONTINUE 

Xi-FREQ-. 5 
X2-FREQ+. 5 
CALL  US TART 
CALL  UBRASE 
CALL  URESET 

CALL  UOAREA (0.00,7.49,0.00,5.71) 

CALL  USET ( "XBOTH") 

CALL  UPSET  ( "XLABEL"  , "FREQUENCY  IN  QHZ 
CALL  USET ( "YBOTH" ) 

CALL  UPSET ( "YLABEL" , "AMPLITUDE; *) 

CALL  USET{ "OWNSCALE") 

CALL  UWINDO (Xl , X2 , 0 . 0 , 2 . 0) 

CALL  UPLOT.l  ( XT ,  DATA  ,4096.0) 

CALL  UPAUSE 
C 

C . . . PLOT  THE  PHASE  SPECTRUM. 


") 


C 


6 


DO  6  1*1,4096 
DATA (I) *0.0 
Yl-AIMAG  v C  ( I ^ ) 

Y2-REAL (C ( I) ) 

Y"Y1**2+Y2**2 
IF(Y.EQ.O)  GOTO  6 
DATA (I) -ATAN2 (Yl,Y2) 

CONTINUE 
CALL  UERASE 

CALL  UDAREA (0.00,7.49,0.00,5.71) 

CALL  UPSET ( "YLABEL" , "PHASE  IN  RADS;") 
CALL  UWINDO (Xl,X2, -4. 0,4.0) 

CALL  UPLOT1 (XT, DATA, 4096.0) 

CALL  UPAUSE 


C 

C... PRINT  REAL  AND  IMAGINARY  PARTS  OF  THE  FREQUENCY 
C... SPECTRUM  IN  FILE  ISPECT. 

C 

PRINT (20, 100)  (XT (I) , REAL (C ( I ) ) ,AIMAG(C(I) ) ,1-1,4097) 

100  FORMAT (1X,2(F7.4,1X,F14.7,1X,F14.7,3X) ) 

125  FORMAT  ( IX ,  F8 . 4 ,  l.X ,  14 ) 

END 

SUBROUTINE  FFT ( X , M , FORWAR) 

Qi.  *******************  **********************  ******************* 

C  FAST  FOURIER  TRANSFORM  VERSION  0.1  30  JULY  1980 

Q*  ******************  ***,f******  ************  ******************** 

C 

C... STANDARD  FAST  FOURIER  TRANSFORM.  X  IS  BOTH  THE  INPUT 
C...AND  THE  OUTPUT  ARRAY  CONTAINING  2**M  COMPLEX  DATA  POINTS. 


If.u 


C. . . FORWAR- . TRUE .  DOES  FORWARD  TRANSFORM,  FORWAR-. FALSE. 
C...DOES  INVERSE  TRANSFORM. 

C 

£************************************************************* 
LOGICAL  FORWAR 
COMPLEX  X,U,W,T 
DIMENSION  X (8200) 

N-2**M 

C 

C...BIT  REVERSAL  SECTION 
C 

NV2-N/2 

NMl-N-1 

J-l 

DO  40  1*1, NMl 

IF  (I.GE.J)  GO  TO  10 

T-X(J) 

X ( J) -X (I) 

X  ( I)  -T 
10  K-NV2 

20  IF  (K.GE.J)  GO  TO  30 

J-J-K 
K-K/2 
GO  TO  20 
30  J-J+K 

40  CONTINUE 

C 

C. . .MULTIPLICATION  SECTION 
C 

PIE*-3. 141592653589 
IF  (.NOT.  FORWAR)  PIE—PIE 
DO  70  L*1,M 
LE«2**L 
LE1-LE/2 

ANGLE-PIE/FLOAT ( LEI ) 

U-(l. 0,0.0) 

W-CMPLX (COS (ANGLE) , SIN (ANGLE) ) 

DO  60  J-l, LEI 

DO  50  I*J, N, LE 
IP-I+LEl 
T*X ( IP) *U 
X(IP)  -X(I)  -T 
X(I)  «X(I)  +T 
50  CONTINUE 

U»U*W 

60  CONTINUE 

?  J  CONTINUE 
C 

C... SCALING  SECTION  -  INVERSE  TRANSFORM  ONLY 
C 

IF  (FORWAR)  RETURN 
SCALE-1. 0/FLOAT(N) 


If.  I 


ooo  noooo  oooo 


DO  80  IJ»1,N 
X (IJ) -X(IJ) * SCALE 
80  CONTINUE 
RETURN 
END 


PROGRAM  INTER ( INPUT, OUTPUT, U9 , UFREQ , TAPE10-U9 , 
+TAPE20-UFREQ, IFREQ, TAPE3 0*1 FREQ) 

...THIS  PROGRAM  INTERPOLATES  BETWEEN  THE  NINE  FREQUENCIES 
.  . . THAT  ARE  THE  OUTPUTS  OF  SHEET1 . 

DIMENSION  TA (4097)  ,A(9)  ,P(9)  ,TP(4097)  ,XT(4097) 

READ (30, 200)  FREQ 

...THE  ARRAY  WHICH  WILL  STORE  THE  INTERPOLATED  FREQUENCY 
...OUTPUT  IS  INITIALIZED  TO  ZERO  SO  THAT  ALL  OF  THE 
...FREQUENCIES  OUTSIDE  OF  THE  PULSE  BANDWIDTH  WILL  BE  ZERO. 

DO  1  1-1,4097 
TA ( I ) -0.0 
TP (I) -0.0 

XT( I) -(FREQ/1024.0) *(I-1) 

1  CONTINUE 

. . . READ  THE  INPUT  DATA  FROM  FILE  U9 . 

DO  2  1-1,9 

READ (10,100)  A ( I) , P ( I) 

2  CONTINUE 
TA ( 7 53 )  -A ( 1) 

TA (821) -A ( 2 ) 

TA (889)  -A{ 3) 

TA( 957) -A ( 4) 

T  A  (1025) -A ( 5 ) 

TA(  1093 ) -A ( 6 ) 

TA ( 1161) -A (7 ) 

TA ( 1229) -A ( 8 ) 

TA (1297) -A ( 9 ) 

TP (753) *P(1) 

TP  (821)  =P(2) 

TP (889) *P ( 3 ) 

TP (957) »P ( 4 ) 


lb? 


TP ( 1025) -P (5) 

TP (1093) -P  (6) 

TP ( 1161) -P  (7) 

TP ( 1229) -P  (8) 

TP ( 1297) *P  (9) 

11-753 

C 

C...THE  LINEAR  INTERPOLATION  BEGINS. 

C 

DO  3  1-753,1297 
12-11+68 
J-I-Il 

TA( I) « ( (TA(l2)-TA(Il) )/(l2-Il) ) *J+TA(I1) 
TP(I)-( (TP (12) -TP ( II) )/ (12-11) ) *J+TP ( II) 

IF ( I . EQ. 12)  11-11+68 
3  CONTINUE 

C 

C... WRITE  THE  INTERPOLATED  SHEETl  OUTPUT  TO  FILE  UFREQ. 
C 

PRINT ( 20 , 125)  (XT (I) ,TA(I) ,TP(I) ,1-1,4097) 

100  FORMAT ( IX , F12 . 8 , 3X , F12 . 8) 

125  FORMAT (IX, 2 (F7 . 4 , IX , F14 . 7 , IX , F14 . 7 , 3X) ) 

200  FORMAT ( IX, F8. 4) 

STOP 

END 


PROGRAM  MULT ( INPUT , OUTPUT , UFREQ , I SPECT , OUTMUL , I FREQ , 
+TAPE1 0 -UFREQ , TAPE2 0 - 1 S  PECT , TAPE3 0 -OUTMUL , TAPE40 - 1 FREQ ) 

C 

C. . . THIS  PROGRAM  TAKES  THE  PULSE  FREQUENCY  SPECTRUM  FROM 
C... PROGRAM  PULSE  AND  THE  INTERPOLATED  FREQUENCY  OUTPUT  FROM 
C... PROGRAM  INTER  AND  MULTIPLIES  THEM  TOGETHER  TO  PRODUCE 
C...THE  OUTPUT  PULSE  FREQUENCY  SPECTRUM. 

C 

DIMENSION  UA ( 4 0 9 7 )  ,UP(4097)  ,FR(4097)  ,FI(4097)  ,XT(4097) 

C 

C . . . READ  THE  INPUT  DATA  FROM  FILES  ISPECT  AND  UFREQ. 

C 

READ (10,125)  (XT ( I) ,UA(I) ,UP(I) ,1-1,4097) 

READ (20 , 125)  (XT(I) ,FR(I) ,FI(I) ,1-1,4097) 

READ (40, 150)  FREQ , TAU 
C 

C... CHANGE  THE  MAGNITUDE  AND  PHASE  PARTS  OF  THE  INTERPOLATED 
C... SPECTRUM  INTO  REAL  AND  IMAGINARY  PARTS. 

C 

DO  5  1*1 , 4097 

IF(UP(I) .GE.2.0)  UP(I) -UP(I) -2.0 
IF (UP ( I) • GE .2,0)  UP(I) »UP(I) -2.0 
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IF  ( UP  ( I)  .GE.  2.0)  UP(I)*UP(I)-2.0 
IF(UP(I) .GE.2.0)  UP (I) «UP(I) -2.0 
5  CONTINUE 

PI*3. 141592654 

DO  1  1*1,4097 
UP (I) -UP (I) *PI 
A-UA  ( I) 

B-UP  ( I) 

CALL  PR (A, B, X, Y) 

UA ( I) *X 
UP (I) *Y 

1  CONTINUE 

...NOW  THE  TWO  FREQUENCY  SPECTRUMS  ARE  MULTIPLIED  TOGETHER. 
DO  2  1*1,4097 

X*  (UA  ( I)  *FR(  I)  )  -  (UP  (I)  *FI  ( I)  ) 

Y* (UP ( I) *FR( I) ) +(UA(I) *FI ( I)  ) 

FR(  I)  -X 
FI ( I) »Y 

2  CONTINUE 

...PLOT  THE  OUTPUT  PULSE  FREQUENCY  SPECTRUM. 

DO  4  1*1 ,4097 
XT  ( I)  M  FREQ/1024 . 0) * ( 1-1) 

UA ( I) *GQRT(FR( I) **2+FI (I) **2) *.009 
4  CONTINUE 

XI *FREQ- . 5 
X~'TREQ4-. 5 
CA,,..  USTART 
CALL  URESET 

TTFPAQP 

CALL  UDAREA(Q.00,7. 49,0.00, 5.71) 

CALL  USET ( "XBOTH” ) 

CALL  USET ( " YBOTH" ) 

CALL  UPSET ("XLABEL", "FREQUENCY  IN  GHZ;") 

CALL  UPSET ( " YLABEL " , "AMPLITUDE ; " ) 

CALL  USET  ( ’’OWNSCALE" ) 

CALL  UWINDO ( XI , X2 , 0 . 0 , 2 . 0 ) 

CALL  UPLOT1 (XT, UA, 4097.0) 

CALL  UPAUSE 

...THE  OUTPUT  PULSE  FREQUENCY  SPECTRUM  IS  O'"  PPUT  TO  FILE 
. . . OUTMUL . 

PRINT (30,100)  (XT ( I ) ,FR(I) ,FI(I) ,1=1,4097) 

100  FORMAT (1X,3(F7.4,1X,F10.3,1X,F9.3,3X) ) 

125  FORMAT (IX, 2 (F7. 4 , IX, F14 . 7 , 1X,F14 .7, 3X) ) 

150  FORMAT (1^ F8.4,1X,I4) 

STOP 

END 
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c 

c************************** ******************************* 

c 

SUBROUTINE  PR(A,P,X,Y) 

C 

C. . .THIS  SUBROUTINE  CONVERTS  "OMPLEX  NUMBERS  IN  THE  POLAR 
C. . .REPRESENTATION  TO  THE  RECTANGULAR  REPRESENTATION. 

C 

PI-3.141592654 
Yl-ABS ( PI/2 . 0-P) 

IF(Y1.LT. .001)  GOTO  1 
Y2-ABS (3.0*PI/2.0-P) 

IF ( Y2 . LT. .001)  GOTO  2 
TP-TAN (P) 

X- ( (A**2) / (1+TP**2) ) **.5 
Y-X*TP 

IF(P.GT. PI/1.0. AND. P.LT. 3. Q*PI/2.0)  GOTO  4 
GOTO  3 

1  X-0.0 
Y-A 
GOTO  3 

2  X-0.0 
Y— A 
GOTO  3 

4  X— X 

Y— Y 

3  RETURN 
END 


PROGRAM  IFFT( INPUT, OUTPUT, OUTMUL, OUTTIM, .APE10 -OUTMUL, 
-f*TAPE20-OUTTlM,  IFREQ , TAPE30-IFREQ) 

C 

C . . . THIS  IROGRAM  TAKES  THE  OUTPUT  FREQUENCY  SPECTRUM  AND 
C... PRODUCES  THE  TIME  DOMAIN  PULSE  VIA  THE  INVERSE  FFT. 

C 

INTEGER  TAU 

DIMENSION  FR(8200) ,FI (8200) 

COMPLEX  C ( 8200 ) 

C 

C . . . READ  THE  INPUT  DATA  FROM  FILE  OUTMUL. 
r 

READ ( 10 , 100 )  (XT, FR(I) ,FI(I) ,1-1,4097) 

READ (30, 150)  FREQ, TAU 
C 
C 

C...THE  FREQUENCY  SPECTRUM  IS  PERIODIC  WITH  A  PERIOD  OF  2  PI. 
C...THE  REAL  PART  OF  THE  SPECTRUM  IS  A  EVEN  FUNCTION,  SO  THE 
C... MIRROR  IMAGE  IS  CREATED  FOR  [PI,2PI].  THE  IMAGINARY  PART 
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...OF  THE  SPECTRUM  IS  AN  ODD  FUNCTION,  SO  THE  NEGATIVE 
...MIRROR  IMAGE  IS  CREATED  FOR  [PI,2PI], 

DO  2  1-1,4095 

FR( 4097  +  1)  »FR( 4097-1) 

FI (4097+1) — FI(4097-I) 

2  CONTINUE 

DO  3  1-1,8192 
C  ( I)  -CMP LX  (FR(  I)  ,  FI  ( I)  ) 

3  CONTINUE 

...CALL  THE  INVERSE  FF J  TO  GET  THE  OUTPUT  PULSE  IN 
...THE  TIME  DOMAIN. 
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CALL  FFT(C,NN, .FALSE. ) 

...PLOT  THE  TIME  FUNCTION 

DO  4  1-1,8192 
FR(  I)  -REAL  (C  ( I)  ) 

FI ( I ) - { 1/ ( FREQ *8 . 0) -4097/ (FREQ* 8 . 0) ) 

4  CONTINUE 
Y-0 .0 

DO  5  1-1,8192 
IF(FR( I) .GT.Y)  Y-FR{ I) 

5  CONTINUE 
IF(Y.LE.l.O)  GOTO  9 

DO  6  1-1,8192 
FR( I ) -FR(I)/Y 

6  CONTINUE 

9  CONTINUE 

DO  15  1-1,4096 
TEMP-FR(I) 

FR(  I ) «FR(8193-I) 

FR(  8193-1) -TEMP 
15  CONTINUE 

Xl-TAU/2+30 
CALL  PLOTS (0. ,0. ,99) 

CALL  USTART 
CALL  UERASE 

CALL  UDAREA{0. 00, 7. 49, 0.00,5.71) 

CALL  USET ( "XBOTH" ) 

CALL  USET  ( " YBOTH" ) 

CALL  UPSET ("XLABEL", "TIME  IN  NANOSEC;") 
CALL  UPSET  (  "  YLABEL"  ,  ".AMPLITUDE ;  "  ) 

CALL  USET( "OWNSCALE") 

CALL  UWI NDO ( - X 1 , X 1 , - 1 . 0 , 1 . 0 ) 

CALL  UPLOT1 (FI, FR, 8192.0) 

CALL  UPAUSE 

CALL  PLOT { 0 . ,0. ,999) 

988  CONTINUE 
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c 

C... OUTPUT  THE  TIME  FUNCTION  TO  FILE  OUTTIM. 

PRINT ( 20 , 125)  (FI ( I) ,FR(I) ,1-1,8192) 

.00  FORMAT ( IX , 3 (F7 . 4 , IX , F10 . 3 , IX , F9 . 3 , 3X) ) 

125  FORMAT ( IX , 4 ( F8 . 1 , 2X , F7 . 4 , 4  X ) ) 

150  FORMAT ( IX ,  F8 . 4 ,  IX ,  1 4 ) 

STOP 

END 

SUBROUTINE  FFT  (X  ,M,  FORWAR) 

C******************************* ********************************* 

C  FAST  FOURIER  TRANSFORM  VERSION  0.1  30  JULY  1980 

C**************************************************************** 

c 

C... STANDARD  FAST  FOURIER  TRANSFORM.  X  IS  BOTH  THE  INPUT 
C...AND  THE  OUTPUT  ARRAY  CONTAINING  2**M  COMPLEX  DATA  POINTS. 

C. . .FORWAR-. TRUE.  DOES  FORWARD  TRANSFORM,  FORWAR-. FALSE.  DOES 
C. . . INVERSE  TRANSFORM. 

C 

C* ****************************************** *********** ********** 

LOGICAL  FORWAR 
COMPLEX  X,U,W, T 
DIMENSION  X ( 8200) 

N-2**M 

C 

C. . .BIT  REVERSAL  SECTION 
C 

NV2-N/2 

NMl-N-1 

J-l 

DO  40  1-1 ,NM1 

IF  (I.GE.J)  GO  TO  10 

T-X(J) 

X(J)  -X(I) 

X(I)  -T 
10  K-NV2 

20  IF  (K.GE.J)  GO  TO  30 

J-J-K 
K-K/2 
GO  TO  20 
30  J-J+K 

40  CONTINUE 

C 

C. . .MULTIPLICATION  SECTION 
C 

PIE— 3.141592653589 
IF  (.NOT.  FORWAR)  PIE— PIE 
DO  70  L-l , M 
LE-2*  *L 
LE1-LE/2 

ANGLE* PIE/FLOAT ( LEI ) 

U»(l. 0,0.0) 


W-CMPLX (COS (ANGLE) , SIN (ANGLE) ) 

DO  60  J-1,LE1 

DO  50  I«J,N,LE 
IP-I+LEl 
T-X(IP)  *U 
X(IP)  -X(I)  -T 
X(I) -X(I) +T 
50  CONTINUE 

U-U*W 

60  CONTINUE 

70  CONTINUE 

C... SCALING  SECTION  -  INVERSE  TRANSPORM  ONLY 

C 

IP  (FORWAR)  RETURN 
SCALE- 1.0 /FLOAT (N) 

DO  80  I J-l , N 
X(IJ)  -X(IJ)  *  SCALE 
80  CONTINUE 

RETURN 
END 


PROGRAM  ENVEL ( INPUT , OUTPUT , OUTTIM , OUTEN , TAPEl 0 -OUTTIM , 
+TAPE20 -OUTEN , OUTM , TAPE3 0 -OUTM , I FREQ , TAPE4  0 • I FREQ ) 

C 

C...THIS  PROGRAM  FINDS  THE  ENVELOPE  OF  A  PULSE  STORED  IN  FILE 
C... OUTTIM.  ALSO  THE  1ST, 2ND  AND  3RD  MOMENTS  ARE  CALCULATED 
C ...  AS  WELL  AS  THE  MEAN  SQUARE  PULSE  WIDTH. 

C 

INTEGER  TAU 

DIMENSION  TO(8192) ,EN(501) ,X(501) 

READ (40,250)  FREQ, TAU 
READ (40, 350)  CORR 
C 

C . . . READ  THE  INPUT  DATA  FROM  FILE  OUTTIM. 

C 

READ (10,100)  (Y,TO(I) ,1-1,8192) 

13-500 

X3-I3 

11-4089 

DO  6  1-4093,4101 
IF(TO(I) .GT.TO(Il) )  Il-I 
6  CONTINUE 

J-Il- (13*4) 

L-I1+ (13*4) 

X2-4097-I1 

C 

C... COMPUTE  THE  PULSE  ENVELOPE. 

C 
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DO  3  1-1,13 
EN(I) -0.0 
3  CONTINUE 

K-l 

DO  1  I-J , L, 8 

X (K) - ( (K-l) - (X3/2) - (X2/8 .0) ) /FREQ 

X (K) -X (K) -CORR 

EN ( K ) - ( TO ( I ) **2) **.5 

K-K+l 

1  CONTINUE 

WRITE (30,325)  I1,J,X(251) 

...PLOT  THE  PULSE  ENVELOPE. 

Xl-TAU/2+30 

CALL  PLOTS ( 0 . ,0. ,99) 

CALL  US TART 
CALL  UERASE 

CALL  UDAR£A(0« 00,7.49, 0.00,5.71) 

CALL  USET("XBOTH") 

CALL  USET( "YBOTH") 

CALL  UPSET ("XLABEL", "TIME  IN  NANOSEC;") 
CALL  UPSET ("YLABEL", "AMPLITUDE;") 

CALL  USET( "OWNSCALE") 

CALL  UWINDO (-Xl , XI, 0.0, 1.0) 

CALL  UPLOTl (X,EN,X3) 

CALL  UPAUSE 

CALL  PLOT (0. ,0. ,999) 

..CALCULATE  THE  FIRST  FOUR  MOMENTS. 

A-0.0 

B-0.0 

C-0.0 

D-0.0 

DO  2  1-1,13 
A- (EN ( I) **2) +A 
B-  (EN  ( I)  **2)  *X(I)  +B 
C- (EN ( I) **2)  * (X ( I) **2)  +C 
D- (EN (I) **2) * (X ( I) **3)+D 

2  CONTINUE 
TMl-B/A 
TM2-C/A 
TM0-A 
TM3-D/A 

RMSP-TM2- (TM1**2) 

SK-TM3- ( 3  *TM2*TMl )  +  ( 2* { TMl **3 ) ) 

SK-0.0 

DO  5  1*1,13 

SK-(EN(I) **2) * ( (X (I) -TMl ) **3) +SK 
5  CONTINUE 

SK-SK/ ( ( RMS  P *  * . 5 )  **3) 


U/t 


o  o  o  onin 


SK-SK/A 


.  ..OTJTPMT  THF  TEMPORAL  MOMENTS  AND  THE  MEAN  SQUARE  PULSE 
...WIDTH  TO  FILE  OUTM. 

WRITE (30/275)  FREQ , TAU 
WRITE (30, 125)  TMO 
WRITE(30,150)  TMl 
WRITE (30,175)  TM2 
WRITE ( 30 , 200)  RMSP 
WRITE (30, 225)  TM3 
WRITE(30, 300)  SK 

...OUTPUT  THE  PULSE  ENVELOPE  TO  FILE  OUTEN. 

WRITE ( 20 , 100)  (X  (I)  ,EN(I)  ,1-1,13) 

100  FORMAT ( IX , 4 (F8 . 1 , 2X , F7 , 4 , 4X) ) 

125  FORMAT (IX, F10 . 4 , 11H  0TH  MOMENT) 

150  FORMAT ( IX , F10 . 4 , 11H  1ST  MOMENT) 

175  FORMAT ( IX, F10 . 4 , 11H  2ND  MOMENT) 

200  FORMAT (1X,F10.4,21H  MEAN  SQ.  PULSE  WIDTH) 

225  FORMAT (IX, FI 5. 4, 11H  3RD  MOMENT) 

250  FORMAT ( IX , F8 , 4 , IX , I 4 ) 

275  FORMAT ( IX , F8 . 4 , 4H  GHZ,3X,I4,8H  NANOSEC) 

300  FORMAT ( IX , FI 0 . 4 , 9H  SKLWNESS) 

325  FORMAT (1X,3H***,3X,I7,3X,I9,3X,F8.2) 

350  FORMAT ( IX, F8. 4) 

STOP 

END 
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The  input  parameters  for  the  Lubble  model  are  stored  in  two 
files. 

I.  File  ELFL 

This  file  contains  the  N(z,x)/Nq(z)  information  of  Figure  2. 

The  N(z,x)/Nq(z)  ratio  remains  constant  for  a  given  x  over  the  entire 
range  of  z.  However,  Nq(z)  changes  according  to  Equation  (30>. 

II.  File  DA2 

This  file  stores  17  input  parameters.  The  input  parameters  are: 

1.  K-576,  the  number  of  staps  in  the  x  direction. 

2.  N-620,  the  number  of  steps  in  the  z  direction. 

3.  H-42.3729,  the  step  size  in  meters  in  the  x  direction. 

4.  TAU-500,  the  step  size  in  meters  in  the  z  direction. 

5.  SIGM-.5,  a  factor  used  in  the  parabolic  equation  solution. 

6.  FREQ-2500,  800,  500,  the  v  ,ve  frequency  in  MHz. 

7.  KSKIP-0,  the  number  of  steps  in  the  z  direction  recorded  on  file 
PASS.  If  KSKIP-0  then  PASS  is  not  read. 

8.  ZM-620,  only  every  ZM  steps  in  the  z  direction  is  stored  on  file 
RESULT . 

9.  KK-1,  only  every  KK  steps  in  the  x  direction  is  recorded  on  file 
RESULT. 

10.  KL-0,  the  number  of  steps  in  the  z  direction  recorded  in  previous 
runs. 

11.  AL-1000,  a  factor  used  in  the  weighting  function. 

12.  BET-20C0,  a  factor  used  in  the  weighting  function. 

13.  Z0-105 .0E+03 ,  a  factor  used  in  the  weighting  function. 


14.  200-195 .0E+03,  a  factor  used  in  the  weighting  function. 

3 

15.  EMAX-2.5E+05,  the  electron  density  in  electrons/cm  at  ZMAX. 
This  parameter  is  named  in  Equation  (30). 

3 

16.  EZ,  not  used,  the  electron  density  in  electrons/cm  at  z-0. 

17.  ZMAX-150.0E+03,  the  height  in  meters  of  the  maximum  in  the 
vertical  electron  density  profile.  This  parameter  is  named 

in  Equation  (30). 

Hq,  the  scale  height,  is  equal  to  110  km  in  Program  SHEET  1. 
The  output  of  SHEET  1  is  stored  in  file  RESULT. 
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.V  • ' fact  Equatorial  ionospheric  scintillation  data  at  1541  5 
Mbs  L-band)  and  3945  5  Mh*  C  band)  >bowmg  uma  shifts  of 
up  to  ona  second  between  similar  fadaa  of  tha  two  signals  ara 
presented  Simpie  modal  computations  show  that  syxtematic  re¬ 
fractive  dTects  due  to  equatorial  pUama  bubble*  tn  tha  partic¬ 
ular  propagation  geometry  may  explain  tha  obaarvad  data. 
Implications  on  equatorial  ionospheric  irregularities  and  scin¬ 
tillation  theory  art  discussed 


Introduction 


Electron  density  irregularities  it.-  tha  nighUima  equatorial 
ionosphere  are  responsible  for  causing  scintillations  of  trans- 
tonusphertcully  propugated  radio  waves  with  frequences  as 
high  as  tha  GHt  range.  Tha  onaat  of  scintillation  has  haan  shown 
to  correspond  to  the  pasaaga  of  targe  depletions  in  ionization 
dansitv  through  the  propagation  path  tYeh  et  al  .  1979;  Ba*u  et 
il  .  1  y HO i  Tha  depletions  have  been  named  plasma  hubbies 
•  McClure  et  al .  1979'  and  ar«  also  asaonatad  with  the  plume 
structures  «een  by  radar  beckscotter  that  extend  vertically 
several  hundred  kilometers  to  tha  topside  of  tha  K  layer  'Wood¬ 
man  and  Lalloi.  1976»  and  the  occurrence  of  "range  spread  F*’ 
■in  lonugrams  obtained  in  tl*a  equatorial  region  'Aarons.  19H2» 

The  scintillations  produced  by  the  irregularities  associated 
with  the  piasma  hubbies  ara  often  quite  severe,  producing  satu¬ 
rated  fluctuations  at  VHK  and  peak-tn-peak  fluctuations  of  up 
to  27  db  at  1541  5  MHs  'L- band » and  H  db  at  J945  3  MHz  'C-handi 
:n  tha  data  analyzed  fur  this  report. 

In  this  paper  we  report  one  aspect  of  the  equatorial  GHz 
scintillation  phenomenon  that  has  not  been  reported  previously 
Scintillation  data  observed  simultaneously  at  L-band  and  C- 
bund  have  been  lound  to  exhibit  time  displacements  of  up  to 
-  <ina  second  between  similar  fades  Previous  workers  have  ob- 
o»rv*j  similar  time  shifts  m  muttifraquenev  observations  of 
radio  star  scintillation  .it  VHK  frequencies  it  low  elevation 
.'•igli-x  'Wild  jnd  Roberts.  1955'  jnd  the  scintillation  of 
pulsars  caused  by  the  tenuous  interstellar  plasma  'Pickett 
.md  Lang.  19?  !•  Shishov  '’971  introduced  a  model  of  the 
interstellar  medium  navmg  two  scales  of  inhomogfneities. 
sma!»  random  scintillation  producing  irregularities  md  large, 
non-r-  'i"ni  ones  which  cause  «m*  refraction,  to  explain  the 
observed  i  quenrv  drift  of  interstellar  scintillation  lent ums 
f.’ole  ir.d  Sice  ! *»H0 »  observed  the  time  shift  in  wideba‘J  'hser- 
vaitons  of  the  qua«:ir  1C  27  3  througn  ihe  tnu-rpl.inetnrv  plasma 
Thoir  interpretation  .■»  that  '.he  dispersion  of  'he  «c»niill  ii«»m 
pattern  s  caused  hv  ddTr.iciion  rrli  K'.inn  from  <h.irp  gr.uiients 
in  the  int»rpl  im-tarv  ••on.iron  densitv  In  is»th  i.ivs  inherent 
refractive  propagation  Hf**«  ts  are  rp^ponsihl#*  h»r  the  oh-a-rvtd 
limp  shift  In  this  r«*t>* in  w«  shall  dti'inpt  •  i  ■  \pl  un  our  ••h«er- 
vMtions  liong  -miil.tr  !im-»  with  onial  ioomiI*  ration  of  he 
geometry  m»l  »h*-  pi-up  igjtiun  path 
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Pita  Collection  and  Processing 


The  data  w**re  collected  m  Ute  -January  and  early  February, 
»**K|  at  A w«* tv, ion  Island  (T  9S  8.  14  4/W.  3CfS  dipt  which  i» 
located  naif  the  ,outh*m  err  it  of  the  Appleton  snomtly  in  F- 
return  mnu.it  ton  Signals  from  MAKJSaT  iIS’Wi  war*  received 
on  I o4 k  ?»  MH*  -Lband)  and  J945  5  MHi  iChandi  The  satellite 
w.ii  viewed  at  un  elevation  of  81*  and  alimuth  of  3&6V  so  the 
propagation  path  wit  a  1  moat  straight  overhead  and  the  pmpo* 
4 niton  v actor  was  within  a  ft w  degrees  of  th#  magnetic  meridian 
puna  The  receiving  antannaa  for  both  L  and  C-band  wan  co¬ 
located  pnma  focus  feed*  illuminating  a  single  3  meter  diameter 
puruboln:  dtah  Tha  signal  to  noma  ratio  aa  dalarminad  from  atrip 
chart  recording*  is  -Itf  <Jb  at  l.-bend  and  -2S  db  at  C-band, 
however,  the  l.  bund  aignal  it  also  "contaminated"  with  un- 
pn  diciubl*  level  change*  which  urt  due  to  changing  traffic  load* 
>n  the  *  ‘  to  l.  lumd  transponder  A  typical  signal  level  "jump"  la 
1  lb  The  analog  AOC  signal  waa  recorded  on  tape  and  later 
digitised  lining  11!  hiu  at  j  rate  of  34  5  Hi  •  029  seconds  per 
>  jn.pl mg  parhHli  and  .lored  on  magnetic  up#  The  receiving  and 
*1*11.1  collect  ion  *\  stem  waa  carefully  checked  to  make  certain  that 
■.he  nim* displacement  between  the  scintillation  pattema  ta  indeed 
due  to  irregular  xiruc'ure  in  the  ionosphere 

The  Urn#  displacement  between  the  diffraction  pattern*  at  the 
two  frequencies  wai  Teaaurvil  by  computing  the  delay  time  r 
it *  the  pea*  nl  the  * row  correlation  between  the  fluctuation*  of 
■he  two  signal*  The  accuracy  of  the  croae*correlation  estimate 
dependent  mi  the  signal  to  noiaa  ratio  on  both  channels,  ana 
•m  i he  ••stent  to  which  ’he  non*  m  correlated  between  the  two 
channel*.  In  addition,  the  unpredictable  signal  level  lumps  in 
:h"  I.  hand  -lata  <  an  cgiw  incorrect  results  whan  tha  peak  to 
P*mk  I  hand  scintillation  level  it  low  ■  «  2  dh  peaa  to  peak)  On 
me  has.*  of  urn  pint  at  *tudie»  we  have  set  a  threshold  scintilla 
Ihim  irnle a  S*  of  IS  at  L-band.  above  which  the  effects  ol  noise  end 
«ign  ii  level  jumpu  are  not  significant  A  farther  constraint  was 
adoed  to  insure  that  the  measured  time  shift*  correspond  to 
-ismior  structure  r,  the  two  signals,  only  shift*  which  roire- 
-pond  *o  j  norm ^ii/i-d  cross-correlation  of  greuter  than  5  have 
consul*  .  cd 
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Result!  and  Discussion 


Wr  have  studied  four  scintillation  event*  which  occurred  in 
th«  time  period  batwetn  2100  and  2400  UTC  on  January  2?  and 
30.  1981  Ascension  Uland  local  tun*  and  UTC  ara  the  urn*.) 
The  time  displacement  h*tw**n  similar  fadaa  at  L-band  and  C- 
hand  is  present  io  varying  dagraar  in  all  of  th*  data  whara  a  rig- 
mftcant  cross-correlation  i>  5)  it  present.  In  Figure  1  we  pre¬ 
sent  example*  of  data  which  show  the  tuna  shift  along  with  th* 
correaponding  croea-correlation  functiona.  Th*  C-band  aignal  ia 
on  top  in  both  caaea.  A  poaitlv*  ahtft  of  th*  croee-correlstlon 
peak  meana  that  th*  C-band  aignal  la  leading  th*  lrb*nd  aignal. 

The  temporal  venation  of  th*  tim*  displacement  is  shown  in 
Figure  2  for  a  34  minute  segment  of  data,  where  each  point  rep- 
r*s«nt*  30  seconds  of  scintillation.  This  data  actually  represents 
two  patches  of  scintillation  activity;  a  brief  initial  patch  of 
intense  scintillation  which  lasted  for  approalmately  sia  minutes 
followed  approximately  2  minutes  later  by  another  patch  which 
lasted  for  over  90  minutes.  Th*  section  of  data  from  2227-2229 
UTC  which  has  been  omitted  corresponds  te  th*  lull  in  scintilla¬ 
tion  activity  Th*  mean  scintillation  tndax  for  th*  entire  seg¬ 
ment  shown  in  Figure  2  is  approximately  1  for  C-band  and  3  for 
L-band.  although  th*  C-band  index  went  as  high  a*  18  and 
L-band  us  high  as  .83.  The  average  peak  cross-correlation  is  .83 
and  th*  average  shift  is  11  seconds  with  an  rtna  shift  of  35 
seconds.  There  appetite  to  be  no  obvious  correlation  between  the 
observed  shifts  and  the  level  of  scintillation  activity. 

On*  of  th*  interesting  characteristics  of  th*  data  became 
evident  when  w*  examined  th*  time  shift  observed  at  th*  "edge*" 
of  th*  patch*!  of  scintillation;  where  w*  ue*  "edge"  to  daecrib* 
the  Unit  and  last  lew  minutas  of  a  patch  of  scintillation.  All  four 
of  (he  patches  w*  studied  had  a  well  daflned  onset  and  deesy  of 
activity  which  is  characterised  by  *  sharp  rise  or  fall  of  scintilla¬ 
tion  index  except  for  th*  second  patch  shown  in  Figure  2  which 
decreased  grad  ally  Th*  time  shift  observed  at  all  four  of  ths 
onsets  is  posit  .*  'C-band  Itsds  Lbandi,  while  th*  shift  observsd 
in  the  last  minutes  of  activity  of  th*  three  event*  which  cessed 
abruptly  is  negative  This  duta  is  summarised  in  Tibi*  1  where 
th*  time  of  onset  and  decay  is  pretanted  along  with  th«  spproti- 
inata  numoer  of  minutes  that  the  time  shift  was  positive  'nega¬ 
tive'  at  the  beginning  -endi  of  an  event. 


Table  1 


DaU* 

Timt  of 
Onset 

#  of  minute*  for 
poHtUvi*  *hift 

Timt  of 
tVciy 

#  ofminut»jfor 
ntqntiv#  nhift 

1  27  H  I 

2127  30 

3  3 

2249  30 

2| 

1  27  8! 

2221.30 

2.5 

2227  00 

3,  „ 

I '  27  81 

22:10.00 

ti’i 

gradual 

:  :)!)■«! 

4 

2254  00 

4 

This  feature  can  lie  seen  in  Figure  2  where  th#  beginning  of 
the  plot  shows  a  positive  shift  which  lasts  lor  -2.5  minutes  and 
'hen  goes  negative  corresponding  to  th*  unset  and  decav  of  the 
first  patch  of  scintillation  Tile  shift  is  again  positive  when  the 
scintillation  resumes  ugain  at  2229 
Since  the  drill  ot  the  nighttime  plasms  is  from  west  to  east, 
presumably,  the  onset  and  decay  of  a  scintillation  patch  rorre 
-pond  to  the  eusi  and  west  edges,  respectively,  of  the  scintillation 
producing  structure  in  'he  ionosphere  In  order  to  espl.un  our 
observ  itions  of  tune  shifts  in  terms  of  what  is  known  about  the 
equatorial  plasma  I'ubh'es  from  radar  and  in  situ  measurements, 
m  tile  next  section  we  present  i  simple  inusiel  ol  refraction  on 
‘.lie  vertlf  i'll  extende*!  , -signs  or  wad.-  "  d  pi  ism  i  hobble* 
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We  model  the  bubble  vail  a*  an  extended  doping  interface 
between  the  baekgound  plasma  and  the  bubble  depletion  where 
the  background  refractive  index  la  n  and  the  refractive  index 
inenle  the  bubble  i»  n  *  An.  Let  u  ray  at  frequency  f  be  incident 
on  the  wall  at  angle  <i  with  reapect  to  the  interface  m  ahown  in 
Fig  -In  We  represent  the  angular  deviation  of  the  ray  from  ita 
initial  direction  by  h  which  is  given  by  elementary  optics* 


<5 


11L 

2  f* 


AN 


cot  a 


(1) 


where  f#  ia  the  background  plasma  frequency  and  AN'N  ii  the 
fractional  change  in  electron  density  serose  the  interface-  The 
direction  of  the  deviation  A  is  such  that,  if  the  ray  is  incident  on  the 
interface  from  a  region  of  higher  refractive  index  (depleted  region! 
the  ray  is  bent  toward  the  intert'aco  while  the  opposite  is  true  for 
incidence  from  a  region  of  lower  refractive  index. 

The  separation  between  two  rays  of  different  frequency  aft  a 
vertical  distance  i  from  the  point  of  refraction  on  the  interface  it 
given  by: 


d  •»  •*  f?  cot  a  (21 

if  we  jx.iome  that  the  bubble  welts  are  field  aligned  with  large 
north-south  extent  and  move  with  a  constant  west  to  east  velocity 
vim.-*),  then  the  time  shift  nsec!  between  similar  fades  will  be 
r  *  *■“  Assuming  a  background  plasma  frequency  of  15  MHx,  a 
height  of  refraction  t-  350  km  and  west  to  east  velocity  v»  100 
nvv  then  for  transmission  si  frequencies  of  1.5  GHx  and  4.0 
(HU  we  have' 


!n 


150  1 


AN. 

N 


I  cot  u 


(3) 


Since  the  propagation  path  is  almost  vertical  and  approximatsly  in 
t  he  magnetic  meridian  plane,  «*  may  be  taken  as  the  approximate 
'till"  ol  the  wall  with  aspect  to  vertical  The  average  magnitude  of 
j  e  is  approximately  3  *K*nda  and  since  the  mean  fractional  change  , 
in  electron  density  must  be  between  0  and  1 .  the  average  tilt  of  the 
will l.i  responsible  for  our  observations  must  be  in  the  range: 

0  ••  a  S  2?  14) 

Recently.  Tsunuda  and  Livingston  1 19811  reported  from  coordl* 
nuted  radar  and  maitu  measurements  that  ail  major  bubbles 
observed  were  at  least  OlFe  depleted  !f  wc  take  this  value  then, 
bused  on  our  observations.  .»n  average  tilt  would  be  u«24*.  The 
maximum  shift  observed  is  about  1  second  which  is  possible  if  the 
refracting  portion  of  the  wall  is  within  d  ;7  of  vertical.  We  now  that 
although  equation  » 1  >  was  derived  for  an  abrupt  interface,  it  is  also 
i  rue  (or  a  gradual  change  in  refractive  index  since  Snell's  Law 
hold*  fur  both  situations. 

So  far  we  have  shown  that  refraction  at  a  bubble  well  can 
produce  time  shifts  of  the  same  magnitude  as  those  observed  in  this 
experiment  Next  we  consider  the  sign  of  the  shift  In  Figure  3b 
and  Jr.  we  show  the  two  possible  configurations  of  the  wait  which 
Will  result  in  a  high  frequency  ray  preceding  a  lower  frequency  ray, 
assuming  a  west  to  east  velocity  Note  t  hat  the  two  configurations 
represent  refraction  on  the  east  wall  of  a  plasma  bubble  A  similar 
picture  can  be  drawn  for  the  case  corresponding  tc  the  west  wuil  of 
the  nubble  and  in  this  case  the  lower  frequency  r-  v  will  precede  the 
high  frequence  r;*v  'for  west  to  east  drift  vein  -,tyi  Thus,  whether 
(be  tiiv  is  incident  on  the  bubble  wall  from  aside  or  outside  of  the 
i *u nine,  as  long  as  (he  ingle  >i  is  not  tor  Urge,  we  observe  the 
fluctuation*  on  the  higher  frequency  first  when  the  rays  are  sepa¬ 
rated  hv  retraction  on  the  east  wall  and  vice  versa  when  the  rays 
are  -wpar.iivd  h*.  refraction  nn  the  west  wall  In  our  data,  the  time 
ihift  al  i  he  verv  beginning  irondufu  patch  of  Humiliation  activity 
"•ai  found  to  persist  from  J-7  minutes.  This  corresponds  to  n  hori¬ 
zontal  displacement  ot  1 2  km  to  42  km  for  tht  scintillation  produc¬ 
ing  patch,  .it.iunnng  a  velocity  of  100  m  »  For  the  simple  Miuatmn 
•lepu  ted  m  Figure  J  with  <•*  «M  thi*  would  i  or  respond  to  a  coher¬ 
ent  vertical  extent  of  fnc  bubble  wall  in  the  rang** 


I7f » 


27  am 


h  •  *!4  km 
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This  idealised  modal  explain*  the  shift  obeerved  at  tha  flnt  and 
laat  faw  minutaa  of  a  acintillation  patch,  however  tha  shift  la  (bund 
to  persist  throughout  tha  patch  which  can  laat  for  over  an  hour.  If 
we  conaidcr  a  typical  "real  Ufa''  bubble  ai  modal  lad  by  Tsunoda  and 
Livingston  1 1981 1  from  radar  and  in  eitu  maaeuremanta,  It  la  dear 
that  tha  structure  oaaoclatad  with  a  tingle  patch  of  acintillation 
activity  can  conaiat  of  multiple  vertically  elongated  wedgatike 
ttructuree  which  attend  from  tha  bottemaida  to  above  tha  peak  of 
tha  F-layer,  with  poeaible  atranger  itructura  beyond  but  near  the 
watt  wall  of  tha  depletion,  Thaae  wadgea  are  often  found  to  tilt  to 
the  weit.  to  a  nearly  vertical  ray  may  paaa  in  and  out  of  the 
depitted  region  eeverat  timea  undergoing  refraction  each  time. 
Thua.  refraction  on  the  boundariee  of  plaama  depletion!  can  alao 
explain  tha  anift  obeerved  in  the  middle  of  a  patch  although  we 
cannot  predict  tha  aign  of  the  dieplacement. 

From  a  allphtly  different  viewpoint,  we  have  atudied  the  ahifta 
obaerved  in  tha  middle  of  tha  acintillation  patch  by  eonatdaring  tha 
simple deterministic  modal  uaadby  Wamikatal.  (1980).  Computa¬ 
tion*  baaed  on  thia  modal  have  revealed  tha  importance  of  ataep. 
vertically  extended  horizontal  gradianta.  which  are  knewn  to  be 
present  maide  the  bubbles  during  their  aerly  development,  in 
producing  enhanced  GHt  scintillation  ieveia.  Wa  have  taken  a 
portion  of  tha  in  situ  profile  shown  in  Figure  2  of  Wemik  ft  al. 

■  19801  which  corresponds  to  tha  aagmant  between  2  and  5  km  from 
tha  bubble  canter,  and  containa  four  dlatinct  irregularities  with 
sharp  gradianta.  Tha  parabolic  equation  ha*  baan  solved  for  a  wava 
propagating  through  a  phase  changing  screen  with  phaee  varia¬ 
tion  proportional  to  the  a lectron  content  fluctuation*  given  by  the 
chosen  segment  of  in  situ  data.  In  Figure  4a  the  ampl  nude  pattern* 
at  a  distance  of  350  km  from  th*  screen  for  L-band  and  C-band 
signal*  respectively.  A  characteristic  pattern  enueed  by  th*  dlf- 
'  fraction  on  th*  sharp  adgas  of  th*  two  outarmoat  irregularities  i  at  f 
approximately  3.75  and  4.6  km  from  tha  cantar  of  tha  bubblel  la 
evident.  Similar  pattama  have  been  obeerved  in  the  Aacenamn 
Island  scintillation  data:  in  aome  caaea.  a  dlatinct  diffraction  pat¬ 
tern  ia  saan  several  time*  in  succession  within  only  a  few  minute*. 

Tha  croM-corrrlation  between  th*  amplitude  at  th*  twe  frequen¬ 
cies  was  computed  separately  for  th*  sect  Ion*  between  2  and  3.6  km 
and  3  6  and  5  25  km  from  tha  cetltar  of  th*  bubble  and  Is  shown  In 
Figure  4b.  In  an*  cee*  th*  croee- correlation  maximum  is  at  *  40  m 
and  in  th*  other  it  is  at  -  20  m.  Thett  dtapiacamenta  correspond  to 
time  shifts  of  -  .4  and  -  2  seconds  if  th*  diffracting  structure 
moves  with  a  velocity  of  100  ms. 


Conclusion* 


In  this  note,  w*  have  presented  equatorial  Ionospheric  scintilla¬ 
tion  data  observed  at  L-band  11541.5  MHs)  and  C-band  1 3945  5 
MH<>  along  a  near  vertical  propagation  path  in  the  magnetic 
meridian  plana  which  showed  time  shifts  of  up  to  on*  ewcolvi 
between  similar  fade*  of  th*  two  signals.  Simple  model  computa¬ 
tions  have  indicated  that  thee*  time  shift*  may  lie  caused  by 
systematic  propagation  effects  such  us  refraction  me  medium  with 
structures  that  ure  consistent  with  th*  vertically  extended  pleama 
bubbles  Indeed,  th*  sign  of  the  observed  tun*  ehlft  and  itsdurmt  ion 
during  the  first  and  th*  lest  few  minutes  of  *  scintillation  patch 
have  further  enabled  us  to  eetimat*  that  the  wtUe  of  the  plasma 
bubble  may  have  coherent  vertical  dimensions  of  a  few  tana  of 
kilometer* 

This  preliminary  study  of  th*  time  shift  aspect  of  the  equatorial 
GHi  scintillation  phenomenon  has  reconfirmed  th*  close  relation¬ 
ship  between  equatorial  plasma  bubble*  and  GHs  acintillation 
Morv  imoort  mily  it  indicate*  that  systematic  propagetion  efTects 
in  addition  (r  random  scattering  play  important  roles  in  producing 
the  observed  scintillating  signal  signature. 
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K'gure  Ii  lb  -top.  middle  i  —  Segments  uf  data  '30  seconda  each 
•  how mg  the  time  milt  between  similar  fades  C  band  te  on  top  m 

l»oth  •« 


f  igure  Ic  bottom  i  -  Normalised  crows  correlation  functions  for 
■he  data  shown  in  Figs  la  lb  A  peak  displaced  to  the  right  meant 
the  t’  Sand  x  urn  o'  *  'riding  trie  L-band  *:gnal. 

Figure  J  —  Temporal  variation  of  the  t..n#  mift  til .  •»,.  cross  corre* 
l.u  urn  pen k  lor  a  14  minute  -wginentofu  u  Each  point  repreaenta 
IU  -econd*  ol  scintillation  data 

Figure  la  •  ielt  •  -  •  <  ienmet  v  of  refraction  at  a  vertically  extended 
gradient  with  incidence  angle  •»  showing  the  deviation  ■>  from  the 
initial  ray  direction 


riwMir*.)b  Jc -center.  right  •  —  Two  configurations  of  an  ionisation 
•vail  which  wi.l  cuu.se  th«  high  frequency  'C-barn;*  ray  to  precede 
•  he  Id*  Irequencv  'L-bundi  ray  when  observed  at  a  single  *uuon 
#n  the  ground  The  v«lu-*itv  's  assumed  to  be  from  W-E. 


figur"  4a  topi  —  The  *:gr.al  amplitude  at  a  distance  of  350  km 
ifirii  i  phase  changing  *.«reen  which  her  phase  fluctuation*  in 
proportion  »o  .-lectron  Jc-rtutv  Huctuattona  given  by  in  situ  data 
<  'Sand  i»  »>n  iop  The  in  situ  data  i»  a  section  of  that  reproduced  in 
Figure  l  ->i  Wt-rnik  **t  ul  '19H0i 

Figure  Ui -bottom*  —  »•*  rreUuon  functions  for  the  llret  and 

ust  tull'oi  th.  dm. i  «hn*r  .n  Fig  5a. 
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1  iSii'i  *  h  impaign  Ilhnoiai.  818U1  •  .J  Austen.  A  \V  WermkVundC  H 

Liu 

i-uu  »uri..i  i » no  spheric  -•vmlil  la*  ion  Jafu  jt  1541  j  \1f(z  <l.dmri<n  and 
■M*  I  MM;  i  ‘  i)  tn*l •  ■•no-.vuig  time  ut  up  In  nr.-i  >*cond  ha*  ween 

irnl'ir  I  idei>>Mh«  ».  •< -u»r.iU  .ire  presented  Minpie  m.Mlt'l  coinputOtian* 
•leiwih.it  r»-i  Ktivi- *-ir*-rtaiiu# '''t'u-.ulnrmi  fu.i-omi  but?Kn'*  in 

the  parti*  ular  ptupag  i  •  gvoruorv  m  ,v  esplum  the  "hserved  *iutj.  lm 
pi i-  at  i»iH  *.n  r*iu.it<a  i  »•  <i  .i  regularities  and  -uni  illat mn  i  neury 

art*  *liscuuHi-*l  .'•(  inti ll.it i*  :  ♦.  ri-lract i«ai.  plasma  bublil*-'*' 
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